基于分數(shù)階導數(shù)模型的巷道圍巖變形失穩(wěn)規(guī)律研究

肖書文

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

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

1 基于分數(shù)階導數(shù)的元件模型

1.1 分數(shù)階導數(shù)定義

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

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

1.2 分數(shù)階導數(shù)蠕變元件模型

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

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

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

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

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

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

2 分數(shù)階導數(shù)蠕變本構模型

2.1 蠕變本構模型

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

圖1 Nishihara蠕變模型

2.2 蠕變本構方程

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

圖2 分數(shù)階Nishihara蠕變模型

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

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

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

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

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

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

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

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

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

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

3 圍巖巖石變形失穩(wěn)規(guī)律分析

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

圖3 三軸流變儀

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

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

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

圖5 本構模型與實驗數(shù)據(jù)的擬合曲線

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

表1 本構方程參數(shù)表

4 本構方程參數(shù)分析

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

5 結語

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

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

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