陳銘，張士誠，柳明，馬新仿，鄒雨時，周彤，李寧，李四海

（1. 中國石油大學（北京），北京 102249；2. 中國石油長城鉆探工程有限公司，北京 100101；3. 中國石化石油勘探開發(fā)研究院，北京 100083）

0 引言

水平井體積改造已成為頁巖氣等致密儲集層油氣開發(fā)的主導技術[1]。該技術常采用段塞加砂、小粒徑支撐劑等工藝，縫內(nèi)支撐劑多為局部支撐，且在主裂縫與分支縫交會處存在“轉角支撐”現(xiàn)象。局部支撐、“轉角支撐”導致支撐劑承受的外應力增大，使支撐劑向巖石縫面的嵌入更顯著，從而引起產(chǎn)能遞減甚至裂縫失效[2-3]。因此，體積改造背景下縫內(nèi)支撐劑的嵌入行為研究及其定量分析具有重要意義。

支撐劑嵌入的研究方法包括實驗研究、數(shù)值模擬和理論模型分析。許多學者針對不同巖石類型開展了支撐劑嵌入實驗研究，并得到關于支撐劑嵌入的基本認識[4-8]。但實驗研究通常針對某種特定類型巖石，研究結果難以推廣。Alramahi等[9]采用有限元方法分析單個支撐劑顆粒與巖石表面的接觸嵌入問題，發(fā)現(xiàn)支撐劑嵌入巖石表面主要發(fā)生于彈塑性變形階段。Deng等[10]采用顆粒流數(shù)值模擬方法研究不同粒徑支撐劑與巖石的相互作用，指出支撐劑濃度越大嵌入深度越小。雖然數(shù)值模擬方法可以求解嵌入過程，但計算量較大，不利于工程設計分析，因此目前應用較為廣泛的仍是支撐劑嵌入理論模型分析方法。

理論模型經(jīng)歷了經(jīng)驗公式、半經(jīng)驗公式和理論解的發(fā)展歷程。Huitt等[11]基于幾何關系建立支撐劑嵌入深度計算公式，該公式中有兩個參數(shù)需要通過實驗研究來擬合確定。Volk等[12]研究支撐劑濃度、大小及分布對嵌入深度的影響，并提出基于實驗結果的經(jīng)驗公式。Li等[13]采用雙球體彈性接觸模型建立了支撐劑嵌入的解析公式，該公式理論基礎完善，但僅考慮了彈性變形，無法分析支撐劑彈塑性嵌入。

本文根據(jù)支撐劑嵌入巖石的力學過程建立支撐劑嵌入巖石的本構方程，該方程考慮支撐劑嵌入的彈性、彈塑性和塑性全過程。結合巖石-支撐劑體系的受力分析，給出支撐劑嵌入深度的計算方法，并將本文模型與彈性模型[13]及 Lacy等[5]、郭建春等[6]、盧聰?shù)萚7]實驗結果進行對比，驗證本文計算方法的可行性和準確性。最后基于本文計算方法，分析支撐劑嵌入深度的影響因素。

1 模型建立

1.1 巖石-支撐劑接觸應力

取支撐劑與巖石的接觸截面A（見圖 1）進行分析，其截面積為At。在該截面上作用有流體壓力和支撐劑接觸應力[14]。根據(jù)受力平衡關系可得：

（1）式中接觸應力表示與巖石接觸的支撐劑承受的應力大小，為分析巖石-支撐劑體系內(nèi)單個支撐劑的受力狀態(tài)，需確定與裂縫面接觸的支撐劑分布狀態(tài)。裂縫內(nèi)支撐劑分布包括 3種形式[15]：單層稀疏分布、單層緊密分布和多層分布（見圖 2）。對于多層分布的支撐劑，受閉合應力影響，支撐劑在縫內(nèi)處于壓實狀態(tài)，因此多層分布狀態(tài)下與巖石接觸的支撐劑可視為單層緊密分布。多層分布的支撐劑接觸應力分析與單層緊密分布支撐劑接觸應力分析相同。

圖1 巖石-支撐劑受力分析圖

圖2 支撐劑分布形式示意圖

采用支撐劑面密度描述支撐劑分布狀態(tài)。支撐劑面密度為單位面積內(nèi)與巖石接觸的支撐劑數(shù)量，即：

受巖石表面粗糙度和支撐劑粒徑分布影響，實際嵌入裂縫表面的支撐劑數(shù)量小于支撐劑總量[16-17]。定義支撐劑嵌入比例ε為支撐劑嵌入數(shù)量與單層支撐劑（與巖石表面接觸的支撐劑）總量之比，即：

單個支撐劑平均嵌入面積為：

聯(lián)立（1）式—（4）式，得到單個支撐劑平均接觸應力為：

（5）式描述了接觸應力與嵌入面積、嵌入量及外載荷的關系。由于外載荷、流體壓力和支撐劑分布為已知量，接觸應力的求解關鍵在于嵌入面積的確定。

1.2 支撐劑嵌入本構方程

支撐劑嵌入巖體的過程可描述為圓球在法向載荷作用下與半無限空間體的接觸力學問題。隨著支撐劑嵌入深度的增加，巖體和支撐劑接觸位置的局部變形經(jīng)歷 3種階段——彈性階段、彈塑性階段和塑性階段（見圖3）。建立支撐劑嵌入過程的本構模型，包括嵌入深度與外載荷的關系、嵌入深度與接觸面積的關系，是分析支撐劑嵌入深度的重要步驟。其中接觸面積為支撐劑嵌入巖石的橫截面積。

圖3 支撐劑嵌入階段示意圖

1.2.1 彈性嵌入

支撐劑與巖石的接觸面積較小時，局部變形滿足彈性變形。根據(jù) Hertz理論[18]，支撐劑嵌入巖石的深度與平均接觸應力的關系為：

廣義彈性模量滿足如下關系式[18]：

彈性嵌入階段，支撐劑嵌入巖石的接觸面積與嵌入深度的關系為[19]：

當支撐劑與巖石的接觸應力達到 1.1倍巖石屈服應力時，巖石開始發(fā)生彈塑性變形[19]。因此，根據(jù)（6）式可得彈性階段最大嵌入深度為：

值得注意的是，根據(jù)（9）式確定的彈性階段最大嵌入深度是理論值。

1.2.2 塑性嵌入

由于彈塑性階段的力學特征參數(shù)值主要以數(shù)值解為主，解析公式需要根據(jù)彈性和塑性解進行插值建立，因此首先描述塑性階段變形，再給出彈塑性嵌入的本構方程。

當接觸應力達到巖石硬度時，巖石開始發(fā)生塑性變形[20]，并維持恒定的接觸應力。材料硬度一般為屈服應力的3倍[21]，因此塑性階段接觸應力為：

塑性階段支撐劑嵌入巖石的接觸面積與嵌入深度的關系為[22]：

巖石發(fā)生塑性變形的臨界嵌入深度的準確確定需借助納米壓痕實驗。為方便理論計算，本文采用固體材料發(fā)生塑性變形的臨界嵌入深度計算方法進行計算。對于固體材料，理論和實驗研究表明，開始發(fā)生塑性變形的臨界嵌入深度一般為彈性階段最大嵌入深度的82.5倍[23-24]。因此，塑性變形的臨界嵌入深度為：

需要注意的是，根據(jù)（12）式確定的發(fā)生塑性變形的臨界嵌入深度是基于均質固體材料假設的理論結果。對于實際巖石，考慮其沉積演化等復雜過程的影響，發(fā)生塑性變形的臨界嵌入深度可能會有變化，準確值需通過納米壓痕實驗確定。

1.2.3 彈塑性嵌入

當dc1≤d≤dc2時，接觸區(qū)域處于彈塑性變形階段[25-27]。彈塑性變形階段屬于純彈性和塑性變形的過渡階段，但力學機理較為復雜，難以得到理論解。采用有限元等數(shù)值方法可以獲得數(shù)值解[21]，但數(shù)值方法計算量較大，不利于工程應用。因此，本文采用插值方法建立彈塑性階段的本構方程，方法簡單，便于計算分析。

根據(jù)Francis[28]的實驗研究結果，彈塑性階段接觸應力可描述為：

根據(jù)[dc1,dc2]端點的應力連續(xù)性條件可以確定系數(shù)a1、a2，由此得到：

由于彈塑性階段為彈性和塑性的過渡階段，因此可假設彈塑性階段接觸面積與嵌入深度的關系為：

采用三次樣條插值法確定c(d)的表達式，令：

連續(xù)性條件為：

求解方程組（17）得到：

將（18）式代入（15）式可得：

為分析插值方法的可靠性，將插值法計算結果與Kogut等[27]的有限元法計算結果對比（見圖 4），可以發(fā)現(xiàn)兩者較為接近，平均誤差為5.71%，說明插值方法具有較高的工程計算精度。

圖4 插值結果和有限元結果對比

1.3 支撐劑嵌入深度求解流程

支撐劑嵌入深度的計算方法取決于支撐劑與巖石的嵌入變形階段，而變形階段也是待確定的，因此采用試算法求解支撐劑嵌入深度。首先根據(jù)巖石和支撐劑基本參數(shù)計算臨界嵌入深度dc1和dc2。假設支撐劑嵌入處于彈性階段，則聯(lián)立（5）式、（6）式和（8）式，得到關于d的非線性方程F1(d)，再采用牛頓法求得彈性嵌入深度，若其值小于等于dc1，則假設正確。否則假設支撐劑嵌入處于彈塑性階段，聯(lián)立（5）式、（14）式和（19）式，得到關于d的非線性方程F2(d)，再采用牛頓法求得彈塑性嵌入深度，若其值小于dc2，則假設正確。否則支撐劑嵌入處于塑性階段，聯(lián)立（5）式、（10）式和（11）式，得到關于d的線性方程F3(d)，從而求得支撐劑嵌入深度。采用Matlab軟件編程實現(xiàn)上述求解過程。

2 模型驗證

為驗證本文模型準確性，將本文模型、彈性模型[13]分別與 Lacy等[5]、郭建春等[6]和盧聰?shù)萚7]實驗結果對比（見圖5）。

彈性模型為Li等[13]推導的支撐劑嵌入深度計算模型，該模型基于雙球體彈性接觸理論解建立。與實驗結果對比發(fā)現(xiàn)，彈性模型計算值小于實驗值，因此引入擬合系數(shù)修正彈性模型。然而，針對不同支撐劑和巖石，需要結合實驗結果確定不同的擬合系數(shù)，且擬合系數(shù)的作用是修正彈性解，會造成對支撐劑嵌入的非彈性力學機制的忽視。而本文模型未引入擬合系數(shù)，避免了這些問題。由圖 5可知，本文模型計算結果與修正后彈性模型及實驗結果較為接近，驗證了本文模型的準確性。因此，本文模型不但更符合支撐劑嵌入的力學機制，而且不需要借助實驗引入修正系數(shù)，計算結果也比較準確，具有可行性和準確性。

圖5 本文模型、彈性模型、修正彈性模型與實驗結果對比

3 支撐劑嵌入深度影響因素

利用本文支撐劑嵌入深度計算方法，以昌吉油田蘆草溝組致密油儲集層 JHW020井參數(shù)為基礎數(shù)據(jù)進行支撐劑嵌入深度影響因素分析。該井巖石和支撐劑基本參數(shù)為：巖石彈性模量為35 210 MPa，巖石泊松比為 0.22，巖石屈服應力為 120 MPa；支撐劑粒徑為0.595/0.297 mm（30/50目），取調(diào)和平均值0.362 mm，支撐劑彈性模量為41 306 MPa，支撐劑泊松比為0.25，支撐劑密度為2.8 kg/m3；閉合應力為60 MPa，生產(chǎn)過程裂縫內(nèi)流體壓力為15 MPa。在對某一影響因素進行分析時，只對該影響因素取不同值，其他參數(shù)均使用基本參數(shù)值。

3.1 支撐劑鋪置濃度

支撐劑鋪置濃度為單位面積的支撐劑質量，根據(jù)支撐劑顆粒密度和體積可得到支撐劑鋪置濃度與面密度的轉化關系：

本文所用支撐劑達到最大面密度時的支撐劑鋪置濃度為0.95 kg/m2。當鋪置濃度大于0.95 kg/m2時，支撐劑呈現(xiàn)多層分布形式，但與巖石表面接觸的支撐劑面密度仍為單層緊密鋪置時的最大面密度，因此繼續(xù)增大濃度不會增加與巖石表面接觸的承壓支撐劑數(shù)量。如圖6所示，當支撐劑濃度大于1.00 kg/m2時，支撐劑嵌入深度和接觸應力不再發(fā)生變化，與前述分析一致。圖6a顯示，不同縫內(nèi)流體壓力條件下，支撐劑濃度0.25 kg/m2時的嵌入深度是1.00 kg/m2時的嵌入深度的2～3倍，即當支撐劑為單層稀疏分布時，支撐劑濃度越大，支撐劑嵌入深度越小。

圖6 不同鋪置濃度的支撐劑嵌入特征

圖6b顯示，在本文算例條件下，不同支撐劑濃度下的接觸應力均大于巖石屈服應力，并小于 3倍巖石屈服應力，表明支撐劑嵌入主要發(fā)生于彈塑性變形階段。同時，流體壓力為0～55 MPa時，接觸應力變化緩慢，流體壓力為55～65 MPa時，接觸應力變化顯著。在閉合應力60 MPa、流體壓力為零的條件下，不同支撐劑濃度下的接觸應力為287～315 MPa，平均值約為閉合應力的 5倍。由于支撐劑嵌入深度較小，接觸面積遠小于閉合應力作用面積，根據(jù)（1）式可以得到接觸應力遠大于閉合應力。需要注意的是，實驗測試支撐劑破碎時的應力為閉合應力，并不是接觸應力，因此接觸應力較高并不能說明支撐劑發(fā)生破碎。

3.2 支撐劑嵌入比例

圖7 不同嵌入比例的支撐劑嵌入特征

圖 7中支撐劑嵌入深度與接觸應力隨流體壓力的變化規(guī)律與圖 6類似。相同條件下，支撐劑嵌入深度和接觸應力隨嵌入比例減小而增大。支撐劑嵌入比例為 1時，承壓支撐劑密度最大，單個顆粒的接觸應力最小，因此嵌入深度最小。嵌入比例受支撐劑粒徑分布影響，粒徑分布分散時嵌入比例較小，而粒徑分布集中的支撐劑嵌入比例會趨近于1，因此采用粒徑均一或分布相對集中的支撐劑有利于減小支撐劑嵌入深度。對于段塞加砂工藝，支撐劑分布會出現(xiàn)局部支撐，嵌入比例會減小，從而增大了嵌入深度。因此，段塞加砂工藝要注重后期支撐劑嵌入的分析，在壓裂設計階段可以提高支撐劑濃度來減小嵌入深度。

3.3 閉合應力與縫內(nèi)流體壓力

當流體壓力大于等于閉合應力時，支撐劑接觸應力為零，不發(fā)生嵌入；隨著流體返排和生產(chǎn)進行，裂縫內(nèi)流體壓力不斷減小，作用在支撐劑上的接觸應力不斷增加，嵌入深度逐漸增大。圖 8顯示，流體壓力為零的條件下，與閉合應力為50 MPa時相比，閉合應力為65 MPa時嵌入深度增大25.6%；在本文算例中閉合應力條件（60 MPa）下，與縫內(nèi)流體壓力為40 MPa時相比，縫內(nèi)流體壓力為零時嵌入深度增大127.5%。因此，高閉合應力情況（如深井壓裂）下支撐劑嵌入會更加顯著；生產(chǎn)過程中流體壓力下降會顯著增大支撐劑嵌入深度，在進行產(chǎn)能遞減分析時不能忽視流體壓力下降導致的支撐劑嵌入問題。

3.4 支撐劑粒徑

圖 9a顯示，相同條件下，減小粒徑會增大支撐劑嵌入深度。圖 9b顯示，不同粒徑的支撐劑主要發(fā)生彈塑性嵌入，接觸應力受支撐劑粒徑影響，粒徑越小，嵌入面積越小，接觸應力則越大。因此，采用小粒徑支撐劑時需要注意支撐劑嵌入引起的縫寬變化。

3.5 巖石與支撐劑彈性模量

圖 10a顯示，減小巖石與支撐劑彈性模量比值會增大支撐劑嵌入深度。對于煤巖、塑性頁巖等，較低的彈性模量勢必帶來較大的嵌入深度，與實驗結果[4-5,29]一致。圖 10b顯示，相同條件下，巖石與支撐劑彈性模量比值越大，接觸應力越大。這是因為相同外載與內(nèi)壓作用下，支撐劑所受的接觸應力與接觸面積成反比，巖石彈性模量較低時支撐劑嵌入深度較大，接觸面積較大，接觸應力減小。

圖8 不同閉合應力的支撐劑嵌入特征

圖9 不同粒徑的支撐劑嵌入特征

圖10 不同巖石與支撐劑彈性模量比值的支撐劑嵌入特征

4 討論

本文支撐劑嵌入深度計算模型可引入到導流能力計算公式中，從而進一步分析導流能力變化。盡管本文對嵌入過程的描述更加接近實際，但本文模型還存在以下主要問題。

①支撐劑嵌入的分析忽略了由于閉合應力過高等因素導致的支撐劑破碎問題。而支撐劑破碎會導致嵌入裂縫面的支撐劑數(shù)量和比例減小，從而加劇支撐劑嵌入[30]。

②在巖石-支撐劑體系的受力分析中，本文引入了支撐劑嵌入比例，僅分析嵌入比例對嵌入深度和接觸應力的影響，未深入探討該參數(shù)的確定方法。實際上，嵌入比例受支撐劑粒徑分布和巖石表面凹凸形貌影響，確定其取值需對巖石表面形態(tài)進行準確刻畫。

5 結論

基于巖石-支撐劑體系受力特征和支撐劑嵌入本構方程提出的支撐劑嵌入深度計算方法考慮了支撐劑彈性—彈塑性—塑性嵌入全過程，與彈性模型相比更加適用于支撐劑嵌入的理論分析和計算，且計算結果較準確。

支撐劑主要發(fā)生彈塑性嵌入。支撐劑單層分布時，鋪置濃度越大，接觸應力越小，嵌入深度越??；多層分布時，與巖石接觸的承壓支撐劑緊密分布，接觸應力和嵌入深度不再隨鋪置濃度增大而變化。采用粒徑均一或分布相對集中的支撐劑有利于減小支撐劑嵌入深度。相同條件下，支撐劑嵌入深度隨閉合應力增大而增大，隨縫內(nèi)流體壓力下降而增大。支撐劑粒徑及巖石與支撐劑彈性模量比值越小，嵌入深度越大。

符號注釋：

a1，a2，c1，c2，c3，c4——待定系數(shù)；A——支撐劑與巖石的接觸截面；Ap——支撐劑與巖石的接觸面積，m2；Ape——單個支撐劑平均嵌入面積，m2；Ap,e，Ap,p，Ap,ep——彈性階段、塑性階段和彈塑性階段支撐劑與巖石的接觸面積，m2；′——彈性階段、塑性階段和彈塑性階段接觸面積對嵌入深度的一階導數(shù)，m；At——截面A總面積，m2；c(d)——待定函數(shù)，1≤c≤2；cp——支撐劑鋪置濃度，kg/m2；d——支撐劑嵌入巖石的深度，m；dc1——彈性階段最大嵌入深度，m；dc2——塑性變形的臨界嵌入深度，m；E′——廣義彈性模量，MPa；Ep，Er——支撐劑、巖石的彈性模量，MPa；F1(d)，F(xiàn)2(d)，F(xiàn)3(d)——關于d的方程；n——單位面積內(nèi)支撐劑嵌入數(shù)量，個/m2；nt——支撐劑面密度，個/m2；N——嵌入巖石的支撐劑總量；Nt——與巖石接觸的支撐劑總數(shù)量；pf——流體壓力，MPa；R——支撐劑顆粒半徑，m；ε——支撐劑嵌入比例，無因次；ρ——支撐劑顆粒密度，kg/m3；σc——閉合應力，MPa；σp——支撐劑接觸應力，MPa；σp,e，σp,p，σp,ep——彈性階段、塑性階段和彈塑性階段的接觸應力，MPa；σy——巖石屈服應力，MPa；υp，υr——支撐劑、巖石的泊松比。

[1]胥云, 陳銘, 吳奇, 等. 水平井體積改造應力干擾計算模型及其應用[J]. 石油勘探與開發(fā), 2016, 43(5): 780-786, 798.XU Yun, CHEN Ming, WU Qi, et al. Stress interference calculation model and its application in volume stimulation of horizontal wells[J]. Petroleum Exploration and Development, 2016, 43(5):780-786, 798.

[2]CHEN D, YE Z H, PAN Z J, et al. A permeability model for the hydraulic fracture filled with proppant packs under combined effect of compaction and embedment[J]. Journal of Petroleum Science and Engineering, 2016, 149: 428-435.

[3]CUI A, GLOVER K, WUST R. Elastic and plastic mechanical properties of liquids-rich unconventional shales and their implications for hydraulic fracturing and proppant embedment: A case study of the Nordberg Member in Alberta Canada[R].Minneapolis: 48th US Rock Mechanics/Geomechanics Symposium,2014.

[4]LACY L, RICKARDS A, ALI S. Embedment and fracture conductivity in soft formations associated with HEC, borate and water-based fracture designs[R]. SPE 38590, 1997.

[5]LACY L, RICKARDS A, BILDEN D. Fracture width and embedment testing in soft reservoir sandstone[J]. SPE Drilling & Completion,1998, 13(1): 25-29.

[6]郭建春, 盧聰, 趙金洲, 等. 支撐劑嵌入程度的實驗研究[J]. 煤炭學報, 2008, 33(6): 661-664.GUO Jianchun, LU Cong, ZHAO Jinzhou, et al. Experimental research on proppant embedment[J]. Journal of China Coal Society,2008, 33(6): 661-664.

[7]盧聰, 郭建春, 王文耀, 等. 支撐劑嵌入及對裂縫導流能力損害的實驗[J]. 天然氣工業(yè), 2008, 28(2): 99-101.LU Cong, GUO Jianchun, WANG Wenyao, et al. Experimental research on proppant embedment and its damage to fractures conductivity[J]. Natural Gas Industry, 2008, 28(2): 99-101.

[8]TAN Y, PAN Z, LIU J, et al. Experimental study of permeability and its anisotropy for shale fracture supported with proppant[J]. Journal of Natural Gas Science & Engineering, 2017, 44: 250-264.

[9]ALRAMAHI B, SUNDBERG M. Proppant embedment and conductivity of hydraulic fractures in shales[R]. ARMA 12-291,2012.

[10]DENG S, LI H, MA G, et al. Simulation of shale-proppant interaction in hydraulic fracturing by the discrete element method[J].International Journal of Rock Mechanics & Mining Sciences, 2014,70(9): 219-228.

[11]HUITT J, MCGLOTHLIN B. The propping of fractures in formations susceptible to propping-sand embedment[R]. API 58-115, 1958.

[12]VOLK L, RAIBLE C, CARROLL H, et al. Embedment of high strength proppant into low-permeability reservoir rock[R]. SPE 9867,1981.

[13]LI K W, GAO Y P, LU Y C, et al. New mathematical models for calculating proppant embedment and fracture conductivity[J]. SPE Journal, 2015, 20(3): 496-507.

[14]ZHANG J J, KAMAMENOV A, ZHU D, et al. Laboratory measurement of hydraulic fracture conductivities in the Barnett shale[R]. SPE 163839-PA, 2014.

[15]CHEN C, MARTYSEVICH V, OCONNELL P, et al. Temporal evolution of the geometrical and transport properties of a fracture/proppant system under increasing effective stress[J]. SPE Journal, 2015, 20(3): 527-535.

[16]MUELLER M, AMRO M. Indentation hardness for improved proppant embedment prediction in shale formations[R]. SPE 174227,2015.

[17]AWOLEKE O, ZHU D, HILL A D. New propped-fracture-conductivity models for tight gas sands[J]. SPE Journal, 2016, 21(5): 1-10.

[18]徐芝綸. 彈性力學簡明教程[M]. 4版. 北京: 高等教育出版社,2013: 228-230.XU Zhilun. Short course of elasticity[M]. 4th ed. Beijing: Higher Education Press, 2013: 228-230.

[19]JOHNSON K L. Contact mechanics[M]. Cambridge: Cambridge University Press, 1985.

[20]BIG-ALABO A, HARRISON P, CARTMELL M P. Contact model for elastoplastic analysis of half-space indentation by a spherical impactor[J]. Computers & Structures, 2015, 151: 20-29.

[21]SONG Z, KOMVOPOULOS K. Elastic-plastic spherical indentation:Deformation regimes, evolution of plasticity, and hardening effect[J].Mechanics of Materials, 2013, 61(8): 91-100.

[22]BARTIER O, HERNOT X, MAUVOISIN G. Theoretical and experimental analysis of contact radius for spherical indentation[J].Mechanics of Materials, 2010, 42(6): 640-656.

[23]STRONGE W J. Impact mechanics[M]. Cambridge: Cambridge University Press, 2000.

[24]STRONGE W J. Contact problems for elasto-plastic impact in multi-body systems[J]. Lecture Notes in Physics, 2000, 551:189-234.

[25]ZHAO Y W, MAIETTA D M, CHANG L. An asperity microcontact model incorporating the transition from elastic deformation to fully plastic flow[J]. Journal of Tribology, 2000, 122(1): 86-93.

[26]BRAKE M R. An analytical elastic-perfectly plastic contact model[J].International Journal of Solids & Structures, 2012, 49(22):3129-3141.

[27]KOGUT L, ETSION I. Elastic-plastic contact analysis of a sphere and a rigid flat[J]. Journal of Applied Mechanics, 2002, 69(5):657-662.

[28]FRANCIS H A. Phenomenological analysis of plastic spherical indentation[J]. Journal of Engineering Materials & Technology Transactions of the ASME, 1976, 98(3): 272-281.

[29]ZHANG J J, OUYANG L, HILL A D, et al. Experimental and numerical studies of reduced fracture conductivity due to proppant embedment in shale reservoirs[R]. SPE 170775, 2014.

[30]HAN J H, WANG J Y, PURI V. A fully coupled geomechanics and fluid flow model for proppant pack failure and fracture conductivity damage analysis[R]. SPE 168617, 2014.