肖書文

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

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

1 基于分?jǐn)?shù)階導(dǎo)數(shù)的元件模型

1.1 分?jǐn)?shù)階導(dǎo)數(shù)定義

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

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

1.2 分?jǐn)?shù)階導(dǎo)數(shù)蠕變元件模型

在元件模型中，常運(yùn)用Newton 黏壺表示黏性變形結(jié)構(gòu)，其本構(gòu)方程如式（2）所示：

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

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

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

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

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

2 分?jǐn)?shù)階導(dǎo)數(shù)蠕變本構(gòu)模型

2.1 蠕變本構(gòu)模型

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

圖1 Nishihara蠕變模型

2.2 蠕變本構(gòu)方程

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

圖2 分?jǐn)?shù)階Nishihara蠕變模型

對圖2 中的模型本構(gòu)方程進(jìn)行分析，其中胡克體的本構(gòu)方程如式（6）所示：

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

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

黏彈性體的本構(gòu)模型如式（7）所示：

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

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

黏塑性體中包含摩擦滑塊，摩擦滑塊主要表示為應(yīng)力閾值點(diǎn)，如式（9）所示：

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

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

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

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

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

圖3 三軸流變儀

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

圖4 不同軸向應(yīng)力加載下蠕變曲線

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

圖5 本構(gòu)模型與實(shí)驗(yàn)數(shù)據(jù)的擬合曲線

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

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

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

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

5 結(jié)語

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

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

3）提出的考慮損傷演化的分?jǐn)?shù)階Nishihara 模型的本構(gòu)方程可以很好地表征巖石的蠕變變形。應(yīng)力、導(dǎo)數(shù)階次、損傷系數(shù)是影響蠕變變形的主要參數(shù)，當(dāng)導(dǎo)數(shù)階次為1時(shí)，分?jǐn)?shù)階模型可退化為Nishihara模型。