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

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

關 鍵 詞：巖石蠕變； 分數階微積分； 蠕變損傷； 黏彈塑性

中圖法分類號： TU45

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2024.11.030

0 引 言

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

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

1 分數階微分理論與模型元件

1.1 Riemann-Liouville型分數階積分

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

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

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

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

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

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

當函數f（t）在0的鄰域內可積時，下式成立：

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）為函數f（t）的Laplace變換。

1.2 分數階黏滯體

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

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

式中：ηγ為Abel黏壺的黏滯系數。

當巖石受恒定應力時，對式（5）兩邊進行分數階積分，可得到Abel黏壺的蠕變方程：

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

對式（6）求導，可得：

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

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

1.3 改進的分數階非線性黏滯體

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

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

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

將式（8）代入式（5）中，非定常Abel黏壺的本構關系可以表示如下：

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

當應力恒定時，對式（9）進行分數階積分，得到非定常Abel黏壺的蠕變方程為

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

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

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

1.4 分數階非線性損傷黏塑性體

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

根據孫鈞^［2］的研究，當加載應力超過巖石材料的損傷應力閾值時，巖石試樣開始發生時效損傷。基于此，根據連續損傷力學理論將損傷變量D定義為^［24］

D=1-σσe（12）

式中：σe為有效應力。

非線性損傷體的應變為

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

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

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

dDdt=Cσ1-DV（14）

式中：D為損傷因子；C、V為與材料性質相關的參數。

假設巖石進入加速蠕變階段的起始時間為損傷開始時間，即ta=0。對上式進行積分，并結合邊界條件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為材料常數。對其進行積分，得到：

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

將式（15）代入式（13）中，得到非線性損傷體的應變為

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

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

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

對于分數階非線性黏塑性體，當σ≥σs時，其本構關系為

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

對式（20）兩邊進行分數階積分，得到分數階非線性黏塑性應變為

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

根據元件組合理論，則分數階非線性損傷黏塑性應變可表示為

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

2 分數階黏彈塑性損傷蠕變模型

2.1 一維蠕變模型

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

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

根據疊加原理，總應變為瞬時彈性應變、黏彈性應變、黏性應變和黏塑性應變4部分之和，即：

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

將式（10）、（15）、（19）代入式（18）中，得到一維應力狀態下的分數階黏彈塑性損傷蠕變本構方程：

ε（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分別為衰減蠕變和加速蠕變階段的初始黏滯系數；η2為等速蠕變階段黏滯系數；γ1、γ2分別為表征衰減蠕變和加速蠕變階段特征的參數，且0lt;γ1lt;1，γ2gt;1。

2.2 三維蠕變模型

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

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

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

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

在三維應力狀態下，總應變εij為

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

根據彈塑性理論，材料內部任意一點的應力張量σij和應變張量εij可表示為

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

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

體積模量和剪切模量滿足以下關系：

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

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

根據彈塑性理論，彈性體的本構關系為

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

則對于彈性體的應變可以表示為

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

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

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

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

黏性體的蠕變方程為

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

根據Perzyna黏塑性理論，三維應力狀態下的黏塑性應變滿足：

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

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

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

h為過應力函數，通常取冪函數形式，即：

hFF0n=FF0n（34）

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

F≥0時，采用相關聯流動法則，即Q=F，于是三維應力狀態下的黏塑性應變的表達式為

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

其中，參數P為

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

一般認為偏應力張量在巖石蠕變中起主要作用，因此引用Drucker-Prager屈服準則，巖石屈服函數為

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

式中：J2為第二偏應力不變量。

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

ε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］凍融循環50次的大理巖單軸蠕變試驗數據對一維蠕變模型的準確性進行初步驗證。大理巖的長期強度為40.72 MPa，蠕變曲線如圖5所示，可以發現：試驗的前兩級應力水平均低于大理巖的長期強度，應變與時間表現為非線性關系，穩態蠕變速率逐漸減小直至趨近于0，采用式（24）的第1式進行描述；當加載應力水平為40.75 MPa和50.94 MPa時，巖石試樣沒有迅速進入加速蠕變階段，蠕變曲線只有瞬時變形階段。衰減和穩定蠕變階段，故仍采用式（24）的第2式進行描述；當應力水平為61.13 MPa時，采用式（24）的第3式進行描述，此時該模型可以全面地描述大理巖典型蠕變曲線的3階段特征。采用非線性最小二乘迭代算法進行擬合，得到模型參數如表1所列。

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

3.2 三維蠕變模型驗證

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

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

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

根據表3中的模型參數，結合式（38）的第3式，繪制鹽巖三軸壓縮蠕變全過程模型曲線和試驗數據對比情況，如圖8所示。由圖可知，模型曲線與試驗結果相關性較好，利用式（38）導出的三維蠕變方程可以很好地描述鹽巖的非線性加速蠕變特征，實現了對巖石蠕變過程3階段的模擬。

4 結 論

為利用巖石蠕變模型準確描述巖體在外荷載作用下的力學響應，本文提出了一種巖石分數階黏彈塑性損傷蠕變模型，并根據相關蠕變試驗數據驗證了該模型的有效性。主要結論如下：

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

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

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

參考文獻：

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

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

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

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

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

［6］ 李娜，于曉要.基于分數階微積分的巖石非定常蠕變本構模型［J］.山東農業大學學報（自然科學版），2021，52（2）：288-292.

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

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

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

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

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

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

［13］關順，王來貴，孫闖.巖石分數階統計損傷流變本構模型［J］.遼寧工程技術大學學報（自然科學版），2021，40（6）：518-524.

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

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

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

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

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

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

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

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

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

［23］陳文玲，趙法鎖，弓虎軍.三軸蠕變試驗中云母石英片巖蠕變參數的研究［J］.巖石力學與工程學報，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］徐衛亞，楊圣奇，褚衛江.巖石非線性黏彈塑性流變模型（河海模型）及其應用［J］.巖石力學與工程學報，2006，25（3）：433-447.

［28］趙越，李磊，閆晗，等.水化-凍融耦合作用下大理巖單軸蠕變力學特性［J］.吉林大學學報（地球科學版），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