基于顆粒接觸理論的巖石壓密階段本構模型

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

(中國礦業大學(北京)應急管理與安全工程學院，北京，100083)

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

1 接觸模型的建立

在沉積學中，根據砂巖內部結構，砂巖主要由碎屑顆粒、雜基、膠結物和孔隙這4個部分組成。其中雜基和膠結物又合稱為填隙物。根據砂巖內部碎屑顆粒數量的不同，砂巖的承載形式主要分為2類：顆粒承載和雜基承載。本文主要研究碎屑顆粒數較多的顆粒承載。CHEN等[19]認為，在成巖過程中，巖石存在壓溶作用，巖石中碎屑顆粒之間呈現凹凸狀的接觸。由于顆粒間凹凸體的存在，巖石在初始加載過程中呈現出非線性特征。

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

本文作者通過分析沉積學中關于砂巖壓密階段的承載規律，提出利用顆粒間不平整接觸面作為本構模型的建模基礎。首先建立能夠反映顆粒間不平整接觸的元件模型，然后利用該元件模型建立巖石壓密階段本構模型。從模型結構來看，本文提出的模型更加簡潔，同時物理意義更加明確。

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

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

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

圖2 壓密階段本構關系的元件模型Fig.2 Element model in compaction stage

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

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

2 G-W模型簡介

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

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

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

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

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

在G-W 模型中，認為粗糙面凸起高度滿足正態分布，其概率密度為

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

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

2.2 粗糙面之間的接觸分析

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

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

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

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

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

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

由于顆粒間存在黏結物，忽略剪應力對接觸面的影響。當一個斜向力作用在接觸面時，假設只考慮法向分量對截面的影響。由于本文重點研究壓密階段巖石的本構關系，故暫未考慮當剪應力較大時接觸面滑移以及巖石破壞的問題。

圖5 轉化坐標后的接觸面示意圖Fig.5 Schematic diagram of contact surface after coordinate transformation

3 壓密階段本構關系

3.1 組合元件的本構關系

當元件兩端產生壓縮量l后，令彈簧體的豎直方向位移為l1，接觸體的豎直方向位移為l2，則有

彈簧中的力P1為

式中：k為彈性系數。根據幾何關系，接觸體的法向位移ln為

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

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

根據串聯結構，令P1=P2，聯立式(10)和(12)，計算得到壓縮量l后組合元件中的力。已知模型尺寸后，便能夠計算得到應力-應變關系。

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

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

對于任意元件，

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

3.2 二分法在計算本構關系中的應用

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

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

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

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

圖6 二分法在計算過程中的應用Fig.6 Application of dichotomy in calculation process

3.3 側向壓力對壓密過程的影響

本文第3.1 節中的受力分析過程僅考慮在單軸壓縮條件下的本構關系，下面將分析側向壓力條件下的本構關系。

首先對元件施加側向壓應力，然后再施加豎直方向的位移。通過側向壓力計算得到接觸面上的側向力Pc，然后計算得到側向壓力在接觸面產生的法向初始位移。在接觸面法向存在預壓力：

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

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

聯立P1=P2，當l已知時即可計算得到單個組合元件的力。在已知模型尺寸后，便能夠計算得到應力-應變關系。

4 模型的驗證

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

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

利用MATLAB數值計算軟件對模型進行編程后，模擬得到模型壓密階段的應力-應變關系。不同砂巖本構關系的實驗結果與數值模擬結果對比如圖7所示。圖7中，R2為模型對實驗結果的擬合優度。

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

本文模型對綠砂巖的擬合優度為0.984，對Berea砂巖的擬合優度為0.991，對Boise砂巖的擬合優度為0.972，對Fontainebleau 砂巖的擬合優度為0.978。本文模型對實驗結果模擬過程中，擬合優度均大于0.95，能夠較準確地描述巖石在壓密階段的非線性變形規律，證明了模型的合理性。

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

表1 模型參數Table1 Model parameters

圖7 不同砂巖本構關系的實驗結果與數值模擬結果對比Fig.7 Comparison of experimental results and numerical simulation results for different sandstone constitutive relations

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

圖8 不同側向壓力條件下壓密階段的應力-應變關系Fig.8 Stress-strain relationship in compaction stage under different lateral pressure conditions

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

5 模型參數分析

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

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

取不同粗糙度ω在無側向壓力的條件下進行數值計算，計算結果如圖10所示。

從圖10可以看出：當粗糙度較小時，隨著應變的增大，應力增長較快；當粗糙度較大時，應力增長的速率較慢。這與接觸面中凸起的接觸密切相關，當ω較小時，說明發生接觸的凸起數量迅速增加，導致應力快速增長。

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

圖10 不同粗糙度條件下的應力-應變關系Fig.10 Stress-strain relationship under different roughness conditions

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

張開度與巖石的孔隙率緊密相關。通過表1與圖7的計算結果可以看出：當巖石孔隙率較大時，需采用較大的張開度才能得到較好的數值計算結果；當巖石孔隙率較小時，需采用較小的張開度才能得到較好的數值計算結果。

圖11 不同張開度條件下的應力-應變關系Fig.11 Stress-strain relationship under different opening conditions

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

6 結論

1)本文建立的接觸元件模型在模擬4種巖石彈性階段本構關系過程中，擬合優度均大于0.95，模型的模擬結果與實驗結果一致，說明模型能夠準確地描述巖石彈性階段的本構關系。

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

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

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