基于分數階導數模型的巷道圍巖變形失穩規律研究

肖書文

（山西澤州天泰錦辰煤業有限公司，山西 晉城 048000）

深部巷道圍巖長期處于高地應力的環境中會發生蠕變變形[1-3]，對于蠕變本構模型的建模研究主要包括經驗模型與組合元件模型，其中組合元件模型由于物理意義清晰而得到了廣泛應用，傳統的元件模型包括Kelvin 模型、Burgers 模型、Nishihara 模型等，它們僅能描述前兩個蠕變變形階段，具有較多局限性[3-5]。本文通過分數階導數的建模方法建立組合元件模型，給出了本構模型的解析解。經數據驗證，基于分數階導數的本構模型可以更好地反映蠕變變形特征，尤其是蠕變的第三階段。

1 基于分數階導數的元件模型

1.1 分數階導數定義

分數階微積分相比于整數階微積分的區別在于積分階次由整數變為了分數，分數階微積分具有多種數學定義方式，其中Riemann-Liouville 分數階微積分定義方法的應用最為廣泛。設f在(0,+∞)上連續分布，并且在[0,+∞] 的任何有限子區間可積，當t＞0，Re(γ)＞0時，稱式（1）為函數f(t)的γ階Riemann-Liouville分數階導數：

基于以上分數階微積分的定義，將傳統的蠕變元件模型進行改進，得到分數階導數的蠕變元件模型。

1.2 分數階導數蠕變元件模型

在元件模型中，常運用Newton 黏壺表示黏性變形結構，其本構方程如式（2）所示：

式中，σ為應力，MPa；η為黏性系數；ε為應變，%；t為時間，h；γ為導數的階次。當γ=1 時，此方程為傳統的Newton 黏壺，用于表示黏性變形狀態。當γ=0 時，此方程為彈簧元件，用于表示理想彈性變形狀態。根據物理意義可知，當0＜γ＜1時，此方程表示分數階黏壺（Abel黏壺），用于表示黏彈性變形狀態。

在蠕變過程中，應力σ為常數，根據Riemann- Liouville 分數階微積分算子理論，將式（2）進行變換得到以時間t為自變量的應變方程如式（3）所示：

在巖石的蠕變過程中通常伴隨著損傷過程，隨著損傷的加劇巖石的變形程度也逐漸增加，最終衍生出宏觀裂隙使巖石發生破壞，在工程中常表現為巷道的變形與失穩的破壞過程。因此，考慮將損傷變量引入到黏性系數中，用于表示黏性系數的劣化過程，將損傷變量定義為負指數形式如式（4）所示：

式中，D為損傷變量；α 為損傷劣化系數，用于表示損傷變化程度。

聯立式（2）、式（3）與式（4），保持應力σ為常數，根據Riemann-Liouville 分數階微積分算子理論，得到以時間t為自變量、包含黏性損傷劣化過程的應變方程如式（5）所示：

2 分數階導數蠕變本構模型

2.1 蠕變本構模型

蠕變過程通常分為初始蠕變階段、穩定蠕變階段與快速蠕變階段，組合原件模型的實際意義是通過不同元件模型的組合搭配可以建立本構方程，以期對上述階段進行描述。其中，Nishihara 模型被廣泛用于描述初始蠕變階段與穩定蠕變階段，如圖1 所示。模型中包括描述彈性應變εe的胡克體、描述黏彈性應變εve的黏彈性體、描述黏塑性應變εvp的黏塑性體。

圖1 Nishihara蠕變模型

2.2 蠕變本構方程

在Nishihara 模型中，將Newton 黏壺替換為分數階Abel黏壺得到圖2所示的分數階Nishihara蠕變模型。

圖2 分數階Nishihara蠕變模型

對圖2 中的模型本構方程進行分析，其中胡克體的本構方程如式（6）所示：

式中，E0為胡克體彈簧元件的彈性模量，GPa。

1.4 統計學處理 采用 SPSS 22.0 軟件進行數據分析。呈正態分布的計量資料以 ±s 表示，兩組間比較采用兩獨立樣本 t 檢驗；呈偏態分布的計量資料以中位數（下四分位數，上四分位數）表示，兩組間比較采用 Mann-Whitney U 檢驗；計數資料以例數和百分數表示，組間比較采用 χ2 檢驗。以是否合并糖尿病微血管病變為因變量，年齡、體質量、腰圍、糖尿病病程、血脂、HbA1c 水平、GA水平、AIR 指數、ACR 指數、INSAUC 和 HOMA-IR指數等為自變量進行 logistic 回歸分析。檢驗水準（α）為 0.05。

黏彈性體的本構模型如式（7）所示：

式中，E1為黏彈性體彈簧元件的彈性模量，GPa；為黏彈性體Abel黏壺的粘性系數。

根據分數階微積分的數學理論，為簡化運算，將Riemann-Liouville 分數階導數替換為Caputo 分數階導數。考慮到初始狀態下，即t=0 時，εVe=0，經過Laplace雙重變換得到式（8）所示的黏彈性體分數階本構方程：

黏塑性體中包含摩擦滑塊，摩擦滑塊主要表示為應力閾值點，如式（9）所示：

式中，σp為摩擦滑塊應力，MPa；σs為閾值應力，MPa。當摩擦滑塊的應力大于閾值應力時，摩擦滑塊發生變形。結合考慮黏性系數劣化過程，并經過 Laplace雙重變化最終得到式（10）所示的黏塑性體本構方程：

綜合考慮圖2 所示的組合元件模型結構，最終得到三階段分數階蠕變模型的本構方程如式（11）所示：

式中，α、β為函數參數；z為函數自變量；k為函數階次。

3 圍巖巖石變形失穩規律分析

為驗證式（13）所示的本構方程的合理性，分析巖石的變形失穩規律，采用山西省晉城市錦辰煤業3 采區巷道圍巖巖石進行單軸應力加載實驗。3 采區整體為單斜構造，水文地質條件中等，地壓正常約為5~9 MPa，3 采區開采條件較好，采區儲量為18.04 Mt，預計巷道的投入使用時間為4~5年，因此，將直接面臨圍巖蠕變變形的影響。本次實驗采用三軸流變儀如圖3所示。

圖3 三軸流變儀

實驗的巖石巖性為粉砂巖，儀器加載的軸向應力為5 MPa、6 MPa、7 MPa、8 MPa、9 MPa，最終得到圖4所示的巖石加載蠕變曲線。

圖4 不同軸向應力加載下蠕變曲線

在實驗中考慮到巷道圍巖所受的水平應力即圍壓通常較小，所以僅進行軸向應力的持續加載。根據圖4 可知，在長時間的軸向應力作用下，巖石發生了蠕變變形。在5~9 MPa 的不同應力加載試驗中，均發生了三階段的蠕變變形，隨著應力的逐漸增加，巖石的變形趨勢更加明顯。其中第一階段的初始蠕變階段主要發生在0 h~25 h 之間，此時巖石發生蠕變變形的速率逐漸減小；第二階段的穩定蠕變階段主要發生在25 ~225 h 之間，巖石發生蠕變的速度保持不變，變形速率為0，巖石進入穩定的變形階段；第三階段的快速蠕變階段主要發生在225 h以后，巖石的蠕變變形速度快速增加，并伴隨著出現宏觀裂隙的產生。在實際工程中，主要原因是在應力的長時間作用下，巖石內部的損傷逐漸積累，在達到破壞的閾值點后發生宏觀破壞，此時巷道圍巖即進入了快速失穩破化階段。相應地在所提出的本構方程中，負指數形式的損傷變量可以很好地表征此過程。最終試驗數據與本構模型具有很好的擬合效果，如圖5所示。

圖5 本構模型與實驗數據的擬合曲線

本構模型的擬合力學參數擬合結果如表1所示，僅改變應力水平即可對不同應力水平下的巖石蠕變變形曲線進行擬合。根據數據可知，隨著應力水平的增加，巖石進入第三階段（快速蠕變階段）的時間逐漸縮短。

表1 本構方程參數表

4 本構方程參數分析

根據式（13）所示的本構模型可知，影響蠕變變形的三個主要參數為應力σ、分數階導數的階次γ、損傷系數α。其中應力σ主要影響巖石不同階段的持續時間，當應力越大時，巖石進入不同變形階段的速度就會加快；分數階導數的階次γ主要影響巖石的變形程度，而不影響巖石進入不同變形階段的速度，當階次逐漸降低時，巖石變形幅度越小，當分數階導數階次為1時，式（13）可以退化為Nishihara 模型。損傷系數α 主要影響巖石第三階段的變形幅度，當損傷系數越大時，巖石進入第三階段的損傷累計速度越快，黏性系數的劣化越逐漸明顯。

5 結語

1）基于分數階微積分理論，將傳統Newton 黏壺改進為分數階Abel 黏壺，并引入負指數形式的黏性系數損傷變量，建立了分數階Nishihara模型的本構方程。

2）對砂巖巖石樣品進行單軸應力加載實驗，得到5~9 MPa 應力下的蠕變變形曲線，隨著應力水平的增加，均發生了三階段蠕變變形，巖石進入第三階段快速蠕變階段的時間逐漸縮短。

3）提出的考慮損傷演化的分數階Nishihara 模型的本構方程可以很好地表征巖石的蠕變變形。應力、導數階次、損傷系數是影響蠕變變形的主要參數，當導數階次為1時，分數階模型可退化為Nishihara模型。