巖石分數(shù)階黏彈塑性損傷蠕變模型研究

摘要：為了準(zhǔn)確描述巖石蠕變變形破壞的全過程特征，基于分數(shù)階微積分理論，提出了一種用于描述巖石初始非線性衰減蠕變階段蠕變變形特征的非定常Abel黏壺；進而根據(jù)連續(xù)損傷力學(xué)理論，建立了考慮蠕變損傷的分數(shù)階非線性損傷黏塑性體，以描述加速蠕變階段的蠕變力學(xué)行為。將非定常Abel黏壺和分數(shù)階非線性損傷黏塑性體與胡克體、黏性體進行串聯(lián)，建立了一種新的巖石分數(shù)階黏彈塑性損傷蠕變模型，并結(jié)合廣義胡克定律及Perzyna黏塑性理論推導(dǎo)出三維應(yīng)力狀態(tài)下的蠕變方程。最后，利用相關(guān)蠕變試驗數(shù)據(jù)反演擬合，對比分析試驗數(shù)據(jù)與模型曲線的相關(guān)性。結(jié)果表明：由該模型導(dǎo)出的蠕變方程不僅可以準(zhǔn)確描述巖石在低應(yīng)力水平下穩(wěn)態(tài)蠕變階段的非線特征，而且能夠反映巖石在高應(yīng)力狀態(tài)下的加速蠕變特征，相關(guān)系數(shù)均在0.96以上，實現(xiàn)了對巖石蠕變過程3階段的模擬。

關(guān) 鍵 詞：巖石蠕變； 分數(shù)階微積分； 蠕變損傷； 黏彈塑性

中圖法分類號： TU45

文獻標(biāo)志碼： A

DOI：10.16232/j.cnki.1001-4179.2024.11.030

0 引 言

巖石的蠕變特性與隧道圍巖、洞室、邊坡等巖體工程的長期穩(wěn)定性密切相關(guān)^［1-5］。巖石蠕變本構(gòu)模型作為定量分析巖石長期強度和變形的途徑之一，其以理論曲線不斷逼近實際工程中巖石的真實蠕變曲線，從而達到模擬巖石蠕變變形特性的目的^［6-8］。當(dāng)前，關(guān)于巖石材料的蠕變特性，研究人員已經(jīng)做了大量的研究和調(diào)查。有研究表明，巖石的蠕變特性受瞬態(tài)應(yīng)力和加載歷史的影響，而分數(shù)階微積分理論具有時間記憶性^［9］，可作為模擬巖石類材料蠕變特性的有力工具，被廣泛應(yīng)用于巖石流變力學(xué)領(lǐng)域。例如，周宏偉^［10］、劉峻松^［11］等采用非定常Abel黏壺來替換西原模型中的Newton黏壺元件，建立了一種鹽巖分數(shù)階蠕變本構(gòu)模型；周瑞鶴等^［12］將分數(shù)階黏塑性體與廣義Kelvin體串聯(lián)，建立了粉砂巖的分數(shù)階蠕變本構(gòu)模型；關(guān)順等^［13］在Koeller彈壺的基礎(chǔ)上，采用統(tǒng)計損傷理論構(gòu)建了一種分數(shù)階統(tǒng)計損傷巖石蠕變本構(gòu)模型，并基于試驗數(shù)據(jù)對模型進行了驗證；于懷昌等^［14］采用Koeller彈壺元件替換整數(shù)階Poynting-Thomson模型中的Newton黏壺元件，構(gòu)建了一種新的巖石分數(shù)階非線性Poynting-Thomson模型；劉文博等^［15］將Abel黏壺的黏滯系數(shù)和分數(shù)階階數(shù)定義為與時間、應(yīng)力相關(guān)的函數(shù)，進而建立了變分數(shù)階非線性蠕變本構(gòu)模型，并用砂巖蠕變試驗結(jié)果進行驗證；胡其志等^［16］將基于分數(shù)階微積分理論構(gòu)建的黏滯體與彈性體串聯(lián)，再與由塑性體和牛頓黏壺并聯(lián)組成的非線性黏塑性體串聯(lián)，建立了四元件蠕變模型。

當(dāng)前，巖石蠕變本構(gòu)模型的研究成果較為豐碩^［17-19］。然而，現(xiàn)有的模型很少能夠同時描述巖石蠕變曲線中衰減蠕變和加速蠕變階段的非線性特征。并且，由于蠕變過程的復(fù)雜性，現(xiàn)有的模型難以描述不同試驗條件下的巖石蠕變特性。綜上所述，本文提出了一種新的巖石分數(shù)階黏彈塑性損傷蠕變模型，并合理推導(dǎo)出三維應(yīng)力狀態(tài)下的蠕變方程，對比分析試驗數(shù)據(jù)與理論曲線的相關(guān)性，以證明該模型可以準(zhǔn)確描述巖石各階段的蠕變曲線特征。

1 分數(shù)階微分理論與模型元件

1.1 Riemann-Liouville型分數(shù)階積分

Riemann-Liouville型由于定義形式簡單，且建模中不需要考慮初始條件^［20］，故較多應(yīng)用于巖石類材料的流變力學(xué)分析中。

對于Riemann-Liouville型定義，當(dāng)tgt;0時，函數(shù)f（t）的γ階分數(shù)階積分算子可以表示為^［21］

Iγf（t）=1Γ（γ）∫tt0（t-ξ）^γ-1f（ξ）dξ（1）

函數(shù)f（t）的γ階微分算子如式（2）所示：

Dγf（t）=1Γ（n-γ）dndtn∫tt0（t-ξ）^n-γ-1f（ξ）dξ（2）

式中：I和D分別表示分數(shù)階積分、微分運算符；γ為分數(shù)階階數(shù)；t為時間；ξ為積分變量；Γ（γ）為伽馬函數(shù)，其表達式為Γ（γ）=∫∞0x^γ-1e^-xdx，Γ（1+γ）=γΓ（γ），γ∈R+，Re（γ）gt;0，n-1lt;γ≤n，n=［γ］。

當(dāng)函數(shù)f（t）在0的鄰域內(nèi)可積時，下式成立：

L［Iγf（t）］=sγF（s）（3）

L［Dγf（t）］=sγF（s）-n-1k=0skD^γ-k-1f（t）t=0（4）

式中：F（s）為函數(shù)f（t）的Laplace變換。

1.2 分數(shù)階黏滯體

研究表明，巖石材料的流變性質(zhì)通常處于理想流體和理想固體之間^［21-22］，為了描述這一特性，將牛頓黏壺本構(gòu)關(guān)系中的一階導(dǎo)數(shù)推廣到分數(shù)階，便得到了Abel黏壺的本構(gòu)關(guān)系，定義為

σ（t）=ηγdγε（t）dtγ（5）

式中：ηγ為Abel黏壺的黏滯系數(shù)。

當(dāng)巖石受恒定應(yīng)力時，對式（5）兩邊進行分數(shù)階積分，可得到Abel黏壺的蠕變方程：

ε（t）=σηγtγΓ（γ+1）（6）

對式（6）求導(dǎo)，可得：

ε·（t）=σηγγt^γ-1Γ（γ+1）（7）

由式（7）可知：當(dāng)γ=1時，Abel黏壺的蠕變速率ε·（t）為常數(shù)σ/η，式（6）可描述等速蠕變階段；當(dāng)γgt;1時，ε·（t）gt;0，式（6）可描述加速蠕變階段；當(dāng)0lt;γlt;1時，ε·（t）lt;0，式（6）可描述衰減蠕變階段。由上述分析結(jié)果可知，Abel黏壺可以描述巖石蠕變曲線的衰減蠕變階段和等速蠕變階段，甚至可通過設(shè)置不同的γ值來調(diào)整蠕變曲線的形狀，描述不同軟硬程度、強度的巖石蠕變特性^［8］。

1.3 改進的分數(shù)階非線性黏滯體

巖石蠕變過程中，其蠕變參數(shù)并非恒定值，而是隨時間的推移不斷發(fā)生變化^［23］。由于Abel黏壺的黏滯系數(shù)是恒定的，故難以準(zhǔn)確地描述巖石的非線性蠕變特性。為了解決這一問題，提出一種非定常Abel黏壺，如圖1所示。假定Abel黏壺的黏滯系數(shù)η與蠕變時間^{［15，22］}滿足如下關(guān)系式：

η=ηγ0e^-λt（8）

式中：ηγ0為改進的分數(shù)階非線性Abel黏壺的初始黏滯系數(shù)，MPa·h；λ為時間影響系數(shù)。

將式（8）代入式（5）中，非定常Abel黏壺的本構(gòu)關(guān)系可以表示如下：

σ（t）=ηγ0e^-λtdγε（t）dtγ（9）

當(dāng)應(yīng)力恒定時，對式（9）進行分數(shù)階積分，得到非定常Abel黏壺的蠕變方程為

ε（t）=σηγ0tγE1，γ+1（λt）（10）

式中：Eα，β（z）為Mittag-Leffler函數(shù)，表達式為

Eα，β（z）=∞k=0zkΓ（αk+β）（z，β∈C；R（α）gt;0）（11）

1.4 分數(shù)階非線性損傷黏塑性體

巖石在衰減和穩(wěn)定蠕變階段內(nèi)只產(chǎn)生可恢復(fù)的黏彈性變形，當(dāng)加載應(yīng)力大于長期強度時，進入加速蠕變階段，巖石材料內(nèi)部的微裂紋開始擴展，損傷不斷積累，最終發(fā)生失穩(wěn)破壞^［22］。因此該階段必須考慮巖石材料內(nèi)部的漸近損傷對其力學(xué)參數(shù)的弱化作用。為此，本文提出了一個非線性損傷體，并與上文中構(gòu)建的分數(shù)階非線性黏滯體進行串聯(lián)，再與塑性體并聯(lián)，構(gòu)建一個考慮蠕變參數(shù)隨時間劣化的非線性損傷黏塑性體，如圖2所示。

根據(jù)孫鈞^［2］的研究，當(dāng)加載應(yīng)力超過巖石材料的損傷應(yīng)力閾值時，巖石試樣開始發(fā)生時效損傷?；诖?，根據(jù)連續(xù)損傷力學(xué)理論將損傷變量D定義為^［24］

D=1-σσe（12）

式中：σe為有效應(yīng)力。

非線性損傷體的應(yīng)變?yōu)?/p>

εvp=σ-σeE1=σE1·D1-D（13）

式中：E1為非線性損傷體的彈性模量。

Kachanov^［25］提出的蠕變損傷變量演化方程如式（14）所示，在巖石材料中得到了廣泛的應(yīng)用。

dDdt=Cσ1-DV（14）

式中：D為損傷因子；C、V為與材料性質(zhì)相關(guān)的參數(shù)。

假設(shè)巖石進入加速蠕變階段的起始時間為損傷開始時間，即ta=0。對上式進行積分，并結(jié)合邊界條件t=ta時，D=0；t=tF-ta時，D=1，得到損傷演化方程：

D=1-1-ttF-ta^11+V（15）

式中：tF為巖石蠕變破壞時間，由巖石蠕變試驗確定；ta為加速蠕變的開始時間，可由Krajcinovic等^［26］提出的蠕變損傷率來確定：

dDdt=pσq（1-D）^-q（16）

式中：p、q為材料常數(shù)。對其進行積分，得到：

ta=1p（q+1）σq（17）

將式（15）代入式（13）中，得到非線性損傷體的應(yīng)變?yōu)?/p>

ε1vp（t）=σE11-ttF-ta^-11+V-1（18）

當(dāng)σ≥σs時，摩擦滑塊啟動，非線性損傷黏塑性體的蠕變方程可以表示為

ε1vp（t）=σ-σsE11-ttF-ta^-11+V-1（19）

對于分數(shù)階非線性黏塑性體，當(dāng)σ≥σs時，其本構(gòu)關(guān)系為

ηγ0e^-λtdγε（t）dtγ+σs=σ（20）

對式（20）兩邊進行分數(shù)階積分，得到分數(shù)階非線性黏塑性應(yīng)變?yōu)?/p>

ε2vp（t）=σ-σsηγ0tγ E1，γ+1（λt）（21）

根據(jù)元件組合理論，則分數(shù)階非線性損傷黏塑性應(yīng)變可表示為

εvp=σ-σsηγ0tγ E1，γ+1（λt）+σ-σsE1-ttF-ta^-11+V-1（22）

2 分數(shù)階黏彈塑性損傷蠕變模型

2.1 一維蠕變模型

圖3為不同應(yīng)力水平下巖石的蠕變曲線特征，其中，εe為彈性應(yīng)變，εve為黏彈性應(yīng)變，εv為黏性應(yīng)變，εvp為黏塑性應(yīng)變，ε·為穩(wěn)態(tài)蠕變速率，σs為長期強度。巖石試樣在加載初期產(chǎn)生的瞬時應(yīng)變可用彈性體描述。隨著加載應(yīng)力的增加，巖石進入衰減蠕變階段，其蠕變速率逐漸減小直至趨近于0，此階段的蠕變應(yīng)變用改進的分數(shù)階非線性黏滯體描述；巖石試樣在經(jīng)歷衰減蠕變階段后進入穩(wěn)態(tài)蠕變階段，其蠕變變形隨時間近似呈線性關(guān)系，該階段的蠕變應(yīng)變由黏性元件描述；當(dāng)加載應(yīng)力小于長期強度時，巖石試樣經(jīng)蠕變曲線的前兩個階段后很快進入加速蠕變階段，故該階段的蠕變變形采用分數(shù)階非線性損傷黏塑性體進行描述。

將上文提出的非定常Abel黏壺和非線性損傷黏塑性體與胡克體、黏性體進行串聯(lián)，得到分數(shù)階黏彈塑性損傷蠕變模型，如圖4所示。

根據(jù)疊加原理，總應(yīng)變?yōu)樗矔r彈性應(yīng)變、黏彈性應(yīng)變、黏性應(yīng)變和黏塑性應(yīng)變4部分之和，即：

ε（t）=εe+εve（t）+εv（t）+εvp（t）（23）

將式（10）、（15）、（19）代入式（18）中，得到一維應(yīng)力狀態(tài)下的分數(shù)階黏彈塑性損傷蠕變本構(gòu)方程：

ε（t）=σE0+σηγ10t^γ1E1，γ1+1（λ1t） （σlt;σs，ε·=0）σE0+σηγ10t^γ1E1，γ1+1（λ1t）+ση2t （σlt;σs，ε·gt;0）σE0+σηγ10t^γ1E1，γ1+1（λ1t）+ση2t+σ-σsηγ30t^γ2 E1，γ2+1（λ2t）+ σ-σsE11-ttF-ta^-11+V-1 （σ≥σs）（24）

式中：E0為瞬時彈性模量；ηγ10、ηγ30分別為衰減蠕變和加速蠕變階段的初始黏滯系數(shù)；η2為等速蠕變階段黏滯系數(shù)；γ1、γ2分別為表征衰減蠕變和加速蠕變階段特征的參數(shù)，且0lt;γ1lt;1，γ2gt;1。

2.2 三維蠕變模型

為了描述巖石在三維應(yīng)力狀態(tài)下的蠕變規(guī)律，根據(jù)徐衛(wèi)亞等^［27］研究成果，采用類比法，將一維應(yīng)力狀態(tài)下的巖石分數(shù)階黏彈塑性損傷蠕變模型推廣到三維情況。為了建立三軸應(yīng)力下巖石分數(shù)階黏彈塑性損傷蠕變方程，首先對巖石蠕變特性作如下假設(shè)^［9］：

（1） 巖石蠕變過程中體積變化是彈性的，蠕變特性主要受偏應(yīng)力張量的影響，即由剪切變形引起。

（2） 巖石試樣為各向同性體。

（3） 巖石蠕變過程中，泊松比保持不變。

在三維應(yīng)力狀態(tài)下，總應(yīng)變εij為

εij=εe_ij+ε^veij+εv_ij+ε^vpij（25）

根據(jù)彈塑性理論，材料內(nèi)部任意一點的應(yīng)力張量σij和應(yīng)變張量εij可表示為

σij=Sij+δijσmεij=eij+δijεm（26）

式中：δij為Kronecker張量；σm、εm分別為球應(yīng)力和球應(yīng)變；Sij、eij分別為偏應(yīng)力張量和偏應(yīng)變張量。

體積模量和剪切模量滿足以下關(guān)系：

K=E02（1+μ）G=E03（1-2μ）（27）

式中：K為體積模量；G為剪切模量；μ為巖石材料的泊松比。

根據(jù)彈塑性理論，彈性體的本構(gòu)關(guān)系為

σm=3KεmSij=2Geij（28）

則對于彈性體的應(yīng)變可以表示為

εe_ij=Sij2G0+13Kσmδij（29）

式中：G0為彈性體的黏彈性剪切模量。

根據(jù)文獻［21，27］中所述的方法，用2G和Sij分別替換一維蠕變模型中的E和σ，得到改進的分數(shù)階非線性黏滯體的三維蠕變方程為

ε^veij=Sij2ηγ10t^γ1E1，γ1+1（λ1t）（30）

黏性體的蠕變方程為

εv_ij（t）=Sij2η2t（31）

根據(jù)Perzyna黏塑性理論，三維應(yīng)力狀態(tài)下的黏塑性應(yīng)變滿足：

ε^vpij=P·〈hFF0n〉Qσij（32）

式中：P為等效參數(shù)；〈·〉為屈服條件判別函數(shù)，即：

〈hFF0n〉=0""" Flt;0hFF0n F≥0（33）

h為過應(yīng)力函數(shù)，通常取冪函數(shù)形式，即：

hFF0n=FF0n（34）

式中：F0為屈服函數(shù)初始參考值，其作用是使表達式量綱統(tǒng)一化，且取F0=1；對于巖土材料，通常取冪指數(shù)n=1。

F≥0時，采用相關(guān)聯(lián)流動法則，即Q=F，于是三維應(yīng)力狀態(tài)下的黏塑性應(yīng)變的表達式為

ε^vpij=P·FFσij（35）

其中，參數(shù)P為

P=1ηγ30t^γ2 E1，γ2+1（λ2t）+12G11-ttF-ta^-11+V-1（36）

一般認為偏應(yīng)力張量在巖石蠕變中起主要作用，因此引用Drucker-Prager屈服準(zhǔn)則，巖石屈服函數(shù)為

F=J2-σs3=σ1-σ3-σs3（37）

式中：J2為第二偏應(yīng)力不變量。

在常規(guī)三軸壓縮蠕變試驗條件下，有σ1gt;σ2=σ3，因此三維應(yīng)力狀態(tài)下巖石分數(shù)階黏彈塑性損傷軸向蠕變方程為

ε11（t）=σ1-σ33G0+σ1+2σ39K+σ1-σ33ηγ10t^γ1E1，γ1+1（λ1t） （σ1-σ3lt;σs，ε·11=0）σ1-σ33G0+σ1+2σ39K+σ1-σ33ηγ10t^γ1E1，γ1+1（λ1t）+σ1-σ33η2t （σ1-σ3lt;σs，ε·11gt;0）σ1-σ33G0+σ1+2σ39K+σ1-σ33ηγ10t^γ1E1，γ1+1（λ1t）+σ1-σ33η2t+ σ1-σ3-σs3ηγ30t^γ2E1，γ2+1（λ2t）+16G11-ttF-ta^-11+V-1 （σ1-σ3≥σs） （38）

3 模型驗證

3.1 一維蠕變模型驗證

利用文獻［28］凍融循環(huán)50次的大理巖單軸蠕變試驗數(shù)據(jù)對一維蠕變模型的準(zhǔn)確性進行初步驗證。大理巖的長期強度為40.72 MPa，蠕變曲線如圖5所示，可以發(fā)現(xiàn)：試驗的前兩級應(yīng)力水平均低于大理巖的長期強度，應(yīng)變與時間表現(xiàn)為非線性關(guān)系，穩(wěn)態(tài)蠕變速率逐漸減小直至趨近于0，采用式（24）的第1式進行描述；當(dāng)加載應(yīng)力水平為40.75 MPa和50.94 MPa時，巖石試樣沒有迅速進入加速蠕變階段，蠕變曲線只有瞬時變形階段。衰減和穩(wěn)定蠕變階段，故仍采用式（24）的第2式進行描述；當(dāng)應(yīng)力水平為61.13 MPa時，采用式（24）的第3式進行描述，此時該模型可以全面地描述大理巖典型蠕變曲線的3階段特征。采用非線性最小二乘迭代算法進行擬合，得到模型參數(shù)如表1所列。

根據(jù)表1中的模型參數(shù)，結(jié)合式（24），繪制大理巖模型曲線和試驗數(shù)據(jù)對比情況，如圖6所示。由圖可知，本文所提出的巖石分數(shù)階黏彈塑性損傷蠕變模型可以很好地反映大理巖在不同應(yīng)力水平下的蠕變曲線特征，且擬合曲線和試驗結(jié)果吻合較好，相關(guān)系數(shù)R²均在0.960以上。

3.2 三維蠕變模型驗證

引入文獻［9］中溫度70 ℃、圍壓為5 MPa的鹽巖三軸壓縮蠕變試驗數(shù)據(jù)進一步對該模型的適用性進行驗證。試驗采用分級加載的方式，來研究不同偏應(yīng)力下鹽巖的蠕變性能。通過試驗數(shù)據(jù)擬合反演得到模型參數(shù)，如表2所列。根據(jù)表2中的模型參數(shù)，結(jié)合式（38）的前兩式，繪制圍壓為5 MPa時的鹽巖三軸蠕變模型曲線和試驗數(shù)據(jù)對比情況，如圖7所示。

表2中，鹽巖三軸壓縮蠕變模型曲線和試驗數(shù)據(jù)相關(guān)系數(shù)R²均在0.960以上。另外，由圖7可知，總體來說，模型曲線與鹽巖穩(wěn)態(tài)蠕變試驗數(shù)據(jù)相關(guān)性較好，表明本文所提出的分數(shù)階黏彈塑性損傷蠕變模型能夠完整地描述巖石在不同階段的蠕變變形特征。

為了進一步驗證模型在描述非線性加速蠕變階段時的準(zhǔn)確性，采用文獻［29］中溫度為25 ℃、軸壓為40 MPa和圍壓為20 MPa的鹽巖三軸壓縮蠕變試驗數(shù)據(jù)，對該模型中的參數(shù)進行擬合反演。值得注意的是，文獻中該圍壓下測得鹽巖的長期強度為19 MPa。其中，鹽巖非線性加速蠕變階段的模型參數(shù)如表3所列?？梢园l(fā)現(xiàn)，其加速蠕變階段模型曲線與試驗數(shù)據(jù)相關(guān)性較好，相關(guān)系數(shù)R²=0.992，驗證了該模型的可靠性。

根據(jù)表3中的模型參數(shù)，結(jié)合式（38）的第3式，繪制鹽巖三軸壓縮蠕變?nèi)^程模型曲線和試驗數(shù)據(jù)對比情況，如圖8所示。由圖可知，模型曲線與試驗結(jié)果相關(guān)性較好，利用式（38）導(dǎo)出的三維蠕變方程可以很好地描述鹽巖的非線性加速蠕變特征，實現(xiàn)了對巖石蠕變過程3階段的模擬。

4 結(jié) 論

為利用巖石蠕變模型準(zhǔn)確描述巖體在外荷載作用下的力學(xué)響應(yīng)，本文提出了一種巖石分數(shù)階黏彈塑性損傷蠕變模型，并根據(jù)相關(guān)蠕變試驗數(shù)據(jù)驗證了該模型的有效性。主要結(jié)論如下：

（1） 基于分數(shù)階微積分理論，提出了一種用于描述初始非線性衰減蠕變特征的分數(shù)階非線性黏滯體。

（2） 提出了一種考慮蠕變損傷的非線性損傷黏塑性體，以描述加速蠕變階段的力學(xué)行為。

（3） 將非定常Abel黏壺和分數(shù)階非線性損傷黏塑性體與胡克體、黏性體進行串聯(lián)，提出了一種新的分數(shù)階黏彈塑性損傷蠕變模型，根據(jù)相關(guān)蠕變試驗數(shù)據(jù)驗證了該模型的有效性。由該模型導(dǎo)出的蠕變方程可以很好地描述巖石的非線性蠕變特征，模型曲線和試驗結(jié)果具有較好的一致性，能夠模擬巖石蠕變曲線各階段特征。

參考文獻：

［1］ 蔡美峰，何滿潮，劉東燕.巖石力學(xué)與工程［M］.北京：科學(xué)出版社，2009：1-5.

［2］ 孫鈞.巖石流變及其工程應(yīng)用［M］.北京：中國建筑工業(yè)出版社，1999.

［3］ 李德宏，于明圓，田大鵬，等.蛇紋大理巖蠕變規(guī)律及非線性蠕變模型研究［J］.地下空間與工程學(xué)報，2023，19（2）：420-427.

［4］ 張俊文，霍英昊.深部砂巖分級增量加卸載蠕變特性［J］.煤炭學(xué)報，2021，46（增2）：661-669.

［5］ 許多，吳世勇，張茹，等.錦屏深部大理巖蠕變特性及分數(shù)階蠕變模型［J］.煤炭學(xué)報，2019，44（5）：1456-1464.

［6］ 李娜，于曉要.基于分數(shù)階微積分的巖石非定常蠕變本構(gòu)模型［J］.山東農(nóng)業(yè)大學(xué)學(xué)報（自然科學(xué)版），2021，52（2）：288-292.

［7］ 劉振，楊圣奇，柏正林，等.循環(huán)加卸載下閃長玢巖蠕變特性及損傷本構(gòu)模型［J］.工程科學(xué)學(xué)報，2022，44（1）：143-151.

［8］ 唐皓，王東坡，裴向軍，等.基于分段模擬的巖石蠕變模型［J］.水文地質(zhì)工程地質(zhì)，2017，44（1）：41-47.

［9］ 張勝利，梁衛(wèi)國，肖寧，等.考慮溫度的鹽巖分數(shù)階黏彈塑性蠕變損傷模型［J］.巖石力學(xué)與工程學(xué)報，2022，41（增2）：3198-3209.

［10］周宏偉，王春萍，段志強，等.基于分數(shù)階導(dǎo)數(shù)的鹽巖流變本構(gòu)模型［J］.中國科學(xué)：物理學(xué)力學(xué)天文學(xué)，2012，42（3）：310-318.

［11］劉峻松，黃海峰，黃敏，等.基于分數(shù)階微積分的巖石蠕變損傷本構(gòu)模型［J］.人民長江，2018，49（7）：81-85.

［12］周瑞鶴，程樺，蔡海兵，等.三軸壓縮分級卸荷條件下粉砂巖蠕變特性及蠕變模型［J］.巖石力學(xué)與工程學(xué)報，2022，41（6）：1136-1147.

［13］關(guān)順，王來貴，孫闖.巖石分數(shù)階統(tǒng)計損傷流變本構(gòu)模型［J］.遼寧工程技術(shù)大學(xué)學(xué)報（自然科學(xué)版），2021，40（6）：518-524.

［14］于懷昌，史廣誠，劉漢東，等.基于分數(shù)階微積分的巖石非線性黏彈性應(yīng)力松弛模型研究［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報，2019，27（1）：194-204.

［15］劉文博，張樹光，陳雷.基于非定常分數(shù)階的巖石時效性蠕變模型［J］.采礦與安全工程學(xué)報，2021，38（2）：388-395.

［16］胡其志，王芝超，丁志剛.基于分數(shù)階導(dǎo)數(shù)的巖石蠕變本構(gòu)模型研究［J］.河南理工大學(xué)學(xué)報（自然科學(xué)版），2021，40（6）：163-168.

［17］咸冬冬，楊圣奇，柏正林，等.閃長玢巖蠕變試驗及黏彈-塑性損傷模型研［J］.人民長江，2023，54（12）：225-232，240.

［18］侯迪，許珍，莫至坤，等.干濕循環(huán)條件下巖石蠕變本構(gòu)模型研究及應(yīng)用［J］.人民長江，2023，54（12）：218-224.

［19］虎積元，盛冬發(fā)，陳泰聰，等.巖石三維非線性黏彈塑性損傷蠕變模型研究［J］.河北工程大學(xué)學(xué)報（自然科學(xué)版），2024，41（2）：36-42.

［20］陳文，孫洪廣，李西成.力學(xué)與工程問題的分數(shù)階導(dǎo)數(shù)建模［M］.北京：科學(xué)出版社，2010.

［21］楊俊濤，宋彥琦，馬宏發(fā)，等.考慮硬化和損傷效應(yīng)的鹽巖蠕變本構(gòu)模型研究［J］.巖土力學(xué)，2023，44（10）：2953-2966.

［22］邵珠山，朱意明，陳浩哲，等.深部硬脆性巖石分數(shù)階蠕變損傷模型研究［J］.應(yīng)用力學(xué)學(xué)報，2024，41（4）：853-860.

［23］陳文玲，趙法鎖，弓虎軍.三軸蠕變試驗中云母石英片巖蠕變參數(shù)的研究［J］.巖石力學(xué)與工程學(xué)報，2011，30（增1）：2810-2816.

［24］LYU C，LIU J，REN Y，et al.Study on very long-term creep tests and nonlinear creep-damage constitutive model of salt rock［J］.International Journal of Rock Mechanics and Mining Sciences，2021，146：104873.

［25］KACHANOV M.Effective elastic properties of cracked solids：critical review of some basic concepts［J］.Applied Mechanics Reviews，1992，45（8）：304-335.

［26］KRAJCINOVIC D，SUMARAC D.A mesomechanical model for brittle deformation processes：part Ⅰ［J］.Journal of Applied Medchanics-Transactions of the Asme，1989，56（1）：51-56.

［27］徐衛(wèi)亞，楊圣奇，褚衛(wèi)江.巖石非線性黏彈塑性流變模型（河海模型）及其應(yīng)用［J］.巖石力學(xué)與工程學(xué)報，2006，25（3）：433-447.

［28］趙越，李磊，閆晗，等.水化-凍融耦合作用下大理巖單軸蠕變力學(xué)特性［J］.吉林大學(xué)學(xué)報（地球科學(xué)版），2023，53（4）：1195-1203.

［29］WU F，LIU J，ZOU Q，et al.A triaxial creep model for salt rocks based on variable-order fractional derivative［J］.Mechanics of Time-Dependent Materials，2021，25（1）：101-118.

（編輯：郭甜甜）

Viscoelastic-plastic damage creep model of rock based on fractional order

HU Jiyuan，SHENG Dongfa，CHEN Taicong，LI Ziheng，YU Hongquan

（Civil Engineering Institute，Southwest Forestry University，Kunming 650224，China）

Abstract：

To accurately characterize the entire process of rock creep deformation and damage，based on the theory of fractional-order calculus，a non-constant Abel viscous element was proposed for characterizing the creep deformation of rocks in the initial nonlinear decay creep stage；furthermore，according to the theory of continuum damage mechanics，a fractional-order nonlinear damage visco-plasticity model that takes account of creep damage was developed to characterize the mechanical behavior during the accelerated creep stage.By connecting the non-constant Abel viscous element and the fractional-order nonlinear damage visco-plasticity model in series with the Hookean body and a viscous element，a new fractional-order viscoelastic-plastic damage creep model for rock was established，and the creep equations in three-dimensional stress state were derived by combining the generalized Hooke′s law and Perzyna′s theory of viscoplasticity.Finally，the correlation between the experimental data and the model curves was comparatively analyzed through inverse fitting using relevant creep test data.The results show that the model-derived creep equations can not only accurately describe the nonlinear characteristics of the steady-state creep phase of the rock at low stress levels，but also reflect the accelerated creep characteristics of the rock at high stress levels，with correlation coefficients all exceeding 0.96，thus achieving simulation of the rock′s three-stage creep process.

Key words：

rock creep； fractional order calculus； creep damage； viscoelastic plasticity