基于顆粒接觸理論的巖石壓密階段本構(gòu)模型

馬秋峰，秦躍平，周天白，楊小彬

(中國礦業(yè)大學(xué)(北京)應(yīng)急管理與安全工程學(xué)院，北京，100083)

砂巖是最重要的油氣儲集層，世界上大部分的油氣資源都儲存在砂巖中。此外，砂巖中常常有煤和各種金屬礦床。因此，對于砂巖力學(xué)性質(zhì)的研究具有重要的意義。由于砂巖內(nèi)部結(jié)構(gòu)的特殊性，傳統(tǒng)Hooke模型不能準(zhǔn)確地描述其應(yīng)力-應(yīng)變關(guān)系，其原因如下：砂巖存在壓密階段，巖石表現(xiàn)出非線性特征；砂巖具有塑性硬化和軟化現(xiàn)象，同樣表現(xiàn)出非線性特征。人們對砂巖的力學(xué)性質(zhì)進(jìn)行了大量研究，發(fā)現(xiàn)砂巖的強(qiáng)度受到孔隙率的影響，孔隙率越高，砂巖的強(qiáng)度越低[1-3]。同時(shí)，砂巖屈服面、硬化函數(shù)、塑性變形同樣受到孔隙率的影響[4-8]。BAUD等[6]分析了不同孔隙率下巖石的單軸抗壓強(qiáng)度，提出了影響巖石強(qiáng)度的含孔隙翼型裂紋模型。KULHAWY[7]認(rèn)為隨著孔隙的壓密，彈性模量迅速增大。AL-HARTHI 等[8]分析了孔隙率對巖石彈性模量的影響，并建立了考慮孔隙影響的彈性模量模型。ZHAO 等[9-10]提出了各向異性應(yīng)力條件下多孔巖石的三維模型，在該模型中，一部分孔隙內(nèi)的變形用基于自然應(yīng)變的胡克定律描述，另一部分孔隙內(nèi)的變形用基于工程應(yīng)變的胡克定律描述。李連崇等[11]完善了上述模型，并利用該模型實(shí)現(xiàn)了FLAC3D開發(fā)。王青元等[12]將巖石峰前變形階段分為2個(gè)部分，分別建立了本構(gòu)模型。曹文貴等[13-14]在統(tǒng)計(jì)損傷模型的基礎(chǔ)上，采用宏觀與微觀相結(jié)合的方法，將空隙巖石分為巖石骨架和空隙2個(gè)部分，建立空隙巖石變形分析模型；該模型不僅反映了空隙巖石的應(yīng)變軟化特征，而且還反映了空隙巖石在初始壓密階段的變形非線性特點(diǎn)。近年來，人們從顆粒接觸入手，對砂巖的力學(xué)性質(zhì)進(jìn)行分析。LESZCZYNSKI[15]利用顆粒的接觸黏聚力，模擬顆粒材料的動載過程。NASSAUER 等[16]提出了多面體顆粒接觸計(jì)算方法。FENG[17]提出了非圓球形接觸模型的計(jì)算方法，為離散元計(jì)算提供了接觸模型。袁康等[18]對巖石壓縮荷載作用下內(nèi)部顆粒組分的力學(xué)響應(yīng)進(jìn)行了研究，從細(xì)觀層面揭示了巖石的破裂機(jī)制。目前，有關(guān)砂巖的壓密階段以及砂巖顆粒接觸的研究較多，而通過接觸理論來描述巖石壓密階段的本構(gòu)關(guān)系的研究較少。本文作者以顆粒接觸為研究基礎(chǔ)，認(rèn)為巖石壓密階段中碎屑顆粒之間的孔隙也被壓密，并且?guī)r石壓密階段非線性特征與顆粒間的接觸有關(guān)。通過引入G-W接觸模型，建立接觸元件模型用于表征巖石中碎屑顆粒的接觸，同時(shí)采用彈簧元件表征顆粒本身的線彈性特征，將二者串聯(lián)，用于描述巖石壓密過程中的應(yīng)力-應(yīng)變關(guān)系。利用建立的串聯(lián)模型計(jì)算巖石壓密階段的本構(gòu)關(guān)系，并與實(shí)驗(yàn)結(jié)果進(jìn)行對比以驗(yàn)證模型的合理性。

1 接觸模型的建立

在沉積學(xué)中，根據(jù)砂巖內(nèi)部結(jié)構(gòu)，砂巖主要由碎屑顆粒、雜基、膠結(jié)物和孔隙這4個(gè)部分組成。其中雜基和膠結(jié)物又合稱為填隙物。根據(jù)砂巖內(nèi)部碎屑顆粒數(shù)量的不同，砂巖的承載形式主要分為2類：顆粒承載和雜基承載。本文主要研究碎屑顆粒數(shù)較多的顆粒承載。CHEN等[19]認(rèn)為，在成巖過程中，巖石存在壓溶作用，巖石中碎屑顆粒之間呈現(xiàn)凹凸?fàn)畹慕佑|。由于顆粒間凹凸體的存在，巖石在初始加載過程中呈現(xiàn)出非線性特征。

在描述壓密階段的本構(gòu)關(guān)系時(shí)，傳統(tǒng)本構(gòu)模型[7-10]往往采用孔隙率作為內(nèi)變量來描述巖石的應(yīng)力-應(yīng)變關(guān)系。盡管模型在一定程度上能夠描述壓密過程，然而此類模型只能反映出壓密階段一部分力學(xué)行為，并沒有從巖石壓密過程的機(jī)理上反映出巖石的非線性變形規(guī)律。例如，模型沒有考慮壓密過程中滲透方向等問題。

本文作者通過分析沉積學(xué)中關(guān)于砂巖壓密階段的承載規(guī)律，提出利用顆粒間不平整接觸面作為本構(gòu)模型的建?；A(chǔ)。首先建立能夠反映顆粒間不平整接觸的元件模型，然后利用該元件模型建立巖石壓密階段本構(gòu)模型。從模型結(jié)構(gòu)來看，本文提出的模型更加簡潔，同時(shí)物理意義更加明確。

傳統(tǒng)的元件模型在巖石力學(xué)行為研究中得到了廣泛應(yīng)用。其中，經(jīng)典西原正夫模型(Nishihara model)被廣泛應(yīng)用于描述巖石蠕變特性。西原正夫模型示意圖如圖1所示。該模型主要包括彈簧元件、黏壺元件以及塑性元件。人們在元件模型的基礎(chǔ)上，提出修正原正夫模型或者新的模型[20-23]，但用于描述巖石壓密階段的模型還未見報(bào)道。本文作者基于上述顆粒接觸理論，在G-W 接觸模型的基礎(chǔ)上，提出巖石壓密階段的元件模型，用于描述巖石的本構(gòu)關(guān)系，如圖2所示。

由圖2可知：該元件存在一個(gè)傾斜面，在該傾斜面上存在滿足赫茲接觸理論的凸起，且凸起高度滿足正態(tài)分布。2個(gè)接觸面之間的平均距離為d0(即接觸面之間的張開度，簡稱張開度)。接觸面與垂直方向的夾角為φ(簡稱接觸面傾角)。接觸面法向方向的力-位移關(guān)系可以采用G-W 模型進(jìn)行計(jì)算。同時(shí)假設(shè)該元件寬度和厚度均為1 m。由于顆粒之間接觸面的角度是隨機(jī)的，并沒有固定的方向，因此，接觸面傾角φ用可于描述大量接觸面接觸方向的整體效果。

圖1 西原正夫模型示意圖Fig.1 Diagram of Nishihara model

圖2 壓密階段本構(gòu)關(guān)系的元件模型Fig.2 Element model in compaction stage

接觸元件模型用于描述壓密階段顆粒間的接觸效應(yīng)。此外，假設(shè)顆粒本身滿足線彈性特征。為了表征顆粒本身的變形規(guī)律，在接觸元件上串聯(lián)1個(gè)彈簧元件(見圖2(b))。

由于壓密階段往往應(yīng)力較小，本文假設(shè)碎屑顆粒不發(fā)生塑性變形，即顆粒一直處于彈性狀態(tài)；由于碎屑顆粒之間存在一定的黏聚力，因此，假設(shè)接觸面不發(fā)生相互滑動，僅考慮接觸面的法向作用力。

2 G-W模型簡介

經(jīng)典G-W 模型用于描述顆粒間接觸面的相互作用。G-W 模型中，粗糙表面被描述為一系列凸起的組合，凸起之間通過赫茲定律進(jìn)行計(jì)算。

G-W 模型假設(shè)如下：1)粗糙面剖面凸起的高度服從正態(tài)分布；2)凸起頂端具有相同的曲率；3)凸起變形相互獨(dú)立，互不影響；4)忽略凸起下部的體積變化。

2.1 G-W模型中的粗糙度與概率密度

G-W 模型中的粗糙表面示意圖如圖3所示。圖中，z(x)為在長度lc范圍內(nèi)橫坐標(biāo)為x處的高度。剖面粗糙度均方根ω(簡稱粗糙度)表示粗糙表面頂點(diǎn)高度與粗糙表面平均高度差值的標(biāo)準(zhǔn)差，其表達(dá)式為

G-W 模型中粗糙面的概率密度示意圖如圖4所示。圖中，黑色區(qū)域表示高度大于z的凸起的概率。概率密度函數(shù)?用于描述剖面某點(diǎn)處高度的概率分布。若以粗糙面平均高度為基準(zhǔn)，在高度為(z,z+dz)范圍內(nèi)出現(xiàn)凸起頂點(diǎn)的概率為p(z，z+dz)，則概率密度函數(shù)?(z)為

在G-W 模型中，認(rèn)為粗糙面凸起高度滿足正態(tài)分布，其概率密度為

圖3 粗糙表面示意圖Fig.3 Diagram of surface roughness

圖4 粗糙面概率密度的求解示意圖Fig.4 Schematic map for probability density solution of rough surface

2.2 粗糙面之間的接觸分析

2個(gè)不同粗糙度的平面相互接觸問題可以等效為剛性平面與可變形的粗糙平面相互接觸的問題。等效后平面的粗糙度ωe可用等效前2 個(gè)平面的粗糙度表示：

式中：ω1和ω2分別為2個(gè)接觸面的粗糙度。

將研究對象取為高度為zs的凸起(zs為高于粗糙面平均高度)。當(dāng)剛性平面與粗糙面平均高度平面的距離為d時(shí)：1)若d＞zs，則凸起與剛性平面不接觸；若d＜zs，則凸起與剛性平面相互接觸，兩者相互重疊的量為zs-d。根據(jù)赫茲定律，凸起與剛性平面的接觸力f(zs)為

式中：E為凸起的彈性模量，β為凸起的曲率。與剛性平面接觸的所有凸起的概率prob(zs＞d)為

因此，單位面積上產(chǎn)生的接觸力P(d)為

式中：N為接觸凸起的數(shù)量。為了便于計(jì)算，假設(shè)當(dāng)接觸面位移為0 mm時(shí)，接觸產(chǎn)生的壓力近似為0 N。轉(zhuǎn)化坐標(biāo)后的接觸面示意圖如圖5所示。以接觸面零位移為坐標(biāo)原點(diǎn)，凸起的平均高度距離起始點(diǎn)的距離視為張開度d0。根據(jù)式(7)轉(zhuǎn)化坐標(biāo)后可得位移為ln時(shí)接觸面的法向力P(ln)為

由于顆粒間存在黏結(jié)物，忽略剪應(yīng)力對接觸面的影響。當(dāng)一個(gè)斜向力作用在接觸面時(shí)，假設(shè)只考慮法向分量對截面的影響。由于本文重點(diǎn)研究壓密階段巖石的本構(gòu)關(guān)系，故暫未考慮當(dāng)剪應(yīng)力較大時(shí)接觸面滑移以及巖石破壞的問題。

圖5 轉(zhuǎn)化坐標(biāo)后的接觸面示意圖Fig.5 Schematic diagram of contact surface after coordinate transformation

3 壓密階段本構(gòu)關(guān)系

3.1 組合元件的本構(gòu)關(guān)系

當(dāng)元件兩端產(chǎn)生壓縮量l后，令彈簧體的豎直方向位移為l1，接觸體的豎直方向位移為l2，則有

彈簧中的力P1為

式中：k為彈性系數(shù)。根據(jù)幾何關(guān)系，接觸體的法向位移ln為

接觸體的接觸力P(ln)計(jì)算公式見式(8)。

對接觸力進(jìn)行分解，得到接觸元件在豎直方向的力P2為

根據(jù)串聯(lián)結(jié)構(gòu)，令P1=P2，聯(lián)立式(10)和(12)，計(jì)算得到壓縮量l后組合元件中的力。已知模型尺寸后，便能夠計(jì)算得到應(yīng)力-應(yīng)變關(guān)系。

通過上述關(guān)系，可計(jì)算得到在任意壓縮量下組合單元中的力。為了求解組合單元中的應(yīng)力-應(yīng)變關(guān)系，需要定義每個(gè)元件的寬度與高度。已知巖石顆粒的彈性模量為Es，則彈簧的彈性系數(shù)k為

式中：S為元件所代表的巖石水平方向的截面積，此處為單位面積；h為元件中彈性部分的高度。

對于任意元件，

式中：h0為元件高度，此處取為1 m。

3.2 二分法在計(jì)算本構(gòu)關(guān)系中的應(yīng)用

式(8)是反映接觸元件力-位移關(guān)系的重要表達(dá)式，但無法對其進(jìn)行積分求解。本文采用數(shù)值計(jì)算的方法，計(jì)算得到近似解。由組合元件的力-位移關(guān)系可得：

假設(shè)總壓縮量l已知，求解組合元件的受力。采用二分法進(jìn)行數(shù)值計(jì)算，如圖6所示。

首先將l分為兩等分(見圖6(a))，令l1=l2=l/2，假設(shè)此時(shí)彈簧元件中的力F1大于接觸元件中的力F2，這說明彈簧元件的壓縮量過大，需進(jìn)行下一步修正。將圖中AB段進(jìn)行等分(見圖6(b))，此時(shí)令l1=l/4，l2= 3×l/4。若此時(shí)F1＞F2，則等分左側(cè)AD段，令l1=l/8，l2= 7×l/8；若F1＜F2，則等分DB段，令l1= 3×l/8，l2= 5×l/8；若此時(shí)F1=F2，則說明原假設(shè)l1=l/4，l2= 3×l/4 就是此時(shí)壓縮量分配的計(jì)算結(jié)果。

依照上述方法，不斷進(jìn)行劃分，直到達(dá)到計(jì)算的精度要求，即|F1-F2|＜R0(其中R0為精度判斷參數(shù))。

圖6 二分法在計(jì)算過程中的應(yīng)用Fig.6 Application of dichotomy in calculation process

3.3 側(cè)向壓力對壓密過程的影響

本文第3.1 節(jié)中的受力分析過程僅考慮在單軸壓縮條件下的本構(gòu)關(guān)系，下面將分析側(cè)向壓力條件下的本構(gòu)關(guān)系。

首先對元件施加側(cè)向壓應(yīng)力，然后再施加豎直方向的位移。通過側(cè)向壓力計(jì)算得到接觸面上的側(cè)向力Pc，然后計(jì)算得到側(cè)向壓力在接觸面產(chǎn)生的法向初始位移。在接觸面法向存在預(yù)壓力：

此后，施加豎直方向的位移。假設(shè)接觸元件豎直位移為Δl2，此時(shí)接觸面的法向位移ln為

通過式(8)計(jì)算得到法向接觸力P(ln)，此時(shí)因豎直位移產(chǎn)生的力P2為

聯(lián)立P1=P2，當(dāng)l已知時(shí)即可計(jì)算得到單個(gè)組合元件的力。在已知模型尺寸后，便能夠計(jì)算得到應(yīng)力-應(yīng)變關(guān)系。

4 模型的驗(yàn)證

下面驗(yàn)證模擬孔隙壓密階段巖石本構(gòu)關(guān)系的元件模型的合理性。本文模型的參數(shù)主要包括顆粒的彈性模量E、單位面積凸起的數(shù)量N、凸起的曲率β、接觸面的張開度d0、接觸面粗糙度ω和接觸模型中的傾角φ。

本文選取4 種不同孔隙率砂巖的實(shí)驗(yàn)結(jié)果[24-25]進(jìn)行對比。由于文獻(xiàn)[24-25]中沒有給出具體的模型參數(shù)，本文參考文獻(xiàn)[26-27]中微觀砂巖顆粒的觀察結(jié)果，推測出4種巖石的參數(shù)，如表1所示。由于孔隙率是確定顆粒間的張開度d0的重要參數(shù)，因此也列于表1中。由于對稱性，顆粒之間接觸面傾角φi的范圍為[0°，90°]。對于顆粒接觸面分布較為均勻的巖石而言，接觸模型中的傾角φ可以通過求解傾角φi的數(shù)學(xué)期望得到，φ= 45o。

利用MATLAB數(shù)值計(jì)算軟件對模型進(jìn)行編程后，模擬得到模型壓密階段的應(yīng)力-應(yīng)變關(guān)系。不同砂巖本構(gòu)關(guān)系的實(shí)驗(yàn)結(jié)果與數(shù)值模擬結(jié)果對比如圖7所示。圖7中，R2為模型對實(shí)驗(yàn)結(jié)果的擬合優(yōu)度。

從圖7可以看出：在巖石壓密階段本構(gòu)關(guān)系表現(xiàn)出明顯的非線性特征。在初始加載過程中，巖石應(yīng)力增長較為緩慢，隨著軸向位移的增大，應(yīng)力增長速率逐漸趨于穩(wěn)定。對比圖7(a)和圖7(c)可知：綠砂巖的壓密階段對應(yīng)的應(yīng)變范圍為0～0.1%，Boise砂巖的壓密階段對應(yīng)的應(yīng)變范圍為0～0.23%，這說明對于不同巖石材料，壓密階段對應(yīng)的應(yīng)變范圍不同。

本文模型對綠砂巖的擬合優(yōu)度為0.984，對Berea砂巖的擬合優(yōu)度為0.991，對Boise砂巖的擬合優(yōu)度為0.972，對Fontainebleau 砂巖的擬合優(yōu)度為0.978。本文模型對實(shí)驗(yàn)結(jié)果模擬過程中，擬合優(yōu)度均大于0.95，能夠較準(zhǔn)確地描述巖石在壓密階段的非線性變形規(guī)律，證明了模型的合理性。

利用Fontainebleau砂巖的模型參數(shù)，對不同側(cè)向壓力條件下進(jìn)行壓縮模擬，模擬過程中分別采用0，5，10，20和30 MPa的側(cè)向壓力。

表1 模型參數(shù)Table1 Model parameters

圖7 不同砂巖本構(gòu)關(guān)系的實(shí)驗(yàn)結(jié)果與數(shù)值模擬結(jié)果對比Fig.7 Comparison of experimental results and numerical simulation results for different sandstone constitutive relations

模擬的加載過程如下：首先施加側(cè)向壓力，垂直方向位移清零后施加垂直壓縮位移。圖8所示為位移過程中豎向的應(yīng)力σ1和應(yīng)變ε1的關(guān)系。

圖8 不同側(cè)向壓力條件下壓密階段的應(yīng)力-應(yīng)變關(guān)系Fig.8 Stress-strain relationship in compaction stage under different lateral pressure conditions

從圖8可以看出：在不同側(cè)向壓力的情況下，隨著側(cè)向壓力的增大，初始壓密階段σ1增長的速率越快。當(dāng)側(cè)壓力接近30 MPa時(shí)，應(yīng)力-應(yīng)變曲線趨近于線性。由此可見，數(shù)值模擬結(jié)果反映了巖石壓密階段的變形特點(diǎn)，與實(shí)驗(yàn)規(guī)律基本一致，這說明模型具有合理性。

5 模型參數(shù)分析

5.1 接觸面粗糙度ω對壓密階段的影響

首先對接觸面粗糙度ω進(jìn)行分析。在本模型中，粗糙度ω表示接觸表面凸起高度距離平均高度的離散程度。不同粗糙度的表面如圖9所示。由圖9可見：表面1的離散度更高，因此表面1的粗糙度更大。

取不同粗糙度ω在無側(cè)向壓力的條件下進(jìn)行數(shù)值計(jì)算，計(jì)算結(jié)果如圖10所示。

從圖10可以看出：當(dāng)粗糙度較小時(shí)，隨著應(yīng)變的增大，應(yīng)力增長較快；當(dāng)粗糙度較大時(shí)，應(yīng)力增長的速率較慢。這與接觸面中凸起的接觸密切相關(guān)，當(dāng)ω較小時(shí)，說明發(fā)生接觸的凸起數(shù)量迅速增加，導(dǎo)致應(yīng)力快速增長。

圖9 不同粗糙度的表面Fig.9 Surfaces with different roughness

圖10 不同粗糙度條件下的應(yīng)力-應(yīng)變關(guān)系Fig.10 Stress-strain relationship under different roughness conditions

5.2 接觸面的張開度d0對壓密階段的影響

張開度與巖石的孔隙率緊密相關(guān)。通過表1與圖7的計(jì)算結(jié)果可以看出：當(dāng)巖石孔隙率較大時(shí)，需采用較大的張開度才能得到較好的數(shù)值計(jì)算結(jié)果；當(dāng)巖石孔隙率較小時(shí)，需采用較小的張開度才能得到較好的數(shù)值計(jì)算結(jié)果。

圖11 不同張開度條件下的應(yīng)力-應(yīng)變關(guān)系Fig.11 Stress-strain relationship under different opening conditions

取不同張開度d0在無側(cè)向壓力的條件下進(jìn)行數(shù)值計(jì)算，計(jì)算結(jié)果如圖11所示。從圖11可以看出：當(dāng)張開度較小時(shí)，在經(jīng)歷較小的應(yīng)變后，應(yīng)力迅速增大；當(dāng)張開度較大時(shí)，在經(jīng)歷較大的應(yīng)變后，應(yīng)力才開始逐漸增長。張開度d0綜合反映了壓密過程的長短，當(dāng)張開度較大時(shí)，需要經(jīng)歷較長的壓密過程，而小的張開度只需要經(jīng)歷較少的壓密過程。

6 結(jié)論

1)本文建立的接觸元件模型在模擬4種巖石彈性階段本構(gòu)關(guān)系過程中，擬合優(yōu)度均大于0.95，模型的模擬結(jié)果與實(shí)驗(yàn)結(jié)果一致，說明模型能夠準(zhǔn)確地描述巖石彈性階段的本構(gòu)關(guān)系。

2)本文提出的接觸元件模型能夠反映圍壓對巖石軸向本構(gòu)關(guān)系的影響規(guī)律，證實(shí)了利用接觸模型用于描述壓密階段本構(gòu)關(guān)系的合理性。

3)巖石顆粒間的張開度對巖石壓密階段具有一定影響。顆粒間的張開度越大，壓密過程經(jīng)歷的應(yīng)變越多；相反，張開度越小，巖石顆粒經(jīng)歷較少的應(yīng)變后能夠進(jìn)入線彈性階段。

4)巖石接觸表面的粗糙度對巖石壓密階段具有一定影響。在相同應(yīng)變條件下，粗糙度越大，應(yīng)力增大越緩慢。