基于分數階導數的凍土損傷蠕變模型研究

閆有民, 單靜云, 劉志義

(1.公路建設與養護技術材料及裝備交通運輸行業研發中心, 甘肅 蘭州 730030；2.甘肅路橋建設集團有限公司, 甘肅 蘭州 730030；3. 甘肅萬泰建設工程有限公司, 甘肅 蘭州 730030；4.蘭州理工大學 土木工程學院, 甘肅 蘭州 730050)

0 引言

凍土路基沉降變形隨時間的持續發展一直以來都是路基病害的關鍵問題之一[1-2]。同時，也是公路工程和鐵路工程穩定性的潛在危害，對高寒高海拔地區的工程建設、運營及養護造成了一系列的影響[3-4]。因此，精確模擬凍土蠕變行為是目前高寒高海拔地區路基工程亟待解決的關鍵性問題之一。

迄今為止，研究者們對凍土蠕變進行了大量的試驗研究，重點考慮了溫度、應力和應變率等因素對凍土時間依賴特性的影響[5-9]。其中，結合蠕變過程中應變率的變化，將凍土蠕變分為衰減蠕變、穩定蠕變和加速蠕變階段[8]。此外，Zhou等[9]指出，應力在未抵達屈服應力之前凍土也會產生不可逆的蠕變變形。在理論模型方面，往往采用元件模型來構建凍土的蠕變模型[10-11]，并通過引入指數函數來實現非線性描述。雖然傳統元件模型為模擬凍土蠕變行為提供了基礎，但由于彈性元件的完全彈性和塑性元件的純粹塑性使模型存在一定的局限性。因此，分數階模型的發展很好地改進了對傳統元件模型的限制，并可通過引入變參數效應來實現黏性損傷描述[12]。Liao等[8]使用黏性損傷與分數階導數建立了溫凍土的蠕變模型，并利用試驗結果對分數階蠕變模型的有效性進行了驗證。Zhou等[13]結合黏塑性理論與分數階導數構建了黏性土的蠕變模型。除此之外，分數階微分作為解釋非常規物理現象的一種數學技巧，可解釋土體的非正交特性[14-15]。但是，在上述凍土蠕變模型中，沒有考慮應力抵達屈服應力之前的黏塑性變形。

然而，在以往的凍土蠕變模型中都忽略了應力抵達屈服應力之前的黏塑性變形，且沒有使用分數階導數來刻畫凍土的塑性流動。因此，本研究在已有試驗結果的基礎上對凍土的損傷效應進行了分析。然后，采用分數階導數和過應力理論獲得了一個可以表征凍土黏塑性行為和損傷效應的黏塑性損傷模型，并采用分數階梯度描述了凍土拋物線屈服準則的非正交特性。最后，結合發展的黏塑性損傷模型、胡克彈簧和Abel阻尼器構建了一個凍土的分數階損傷蠕變模型，并使用高溫凍土的蠕變試驗結果驗證了模型的有效性。

1 Caputo分數階導數

分數階微分的定義有多種方法，而在構建物理數學模型中應用較為廣泛的有：Riemann-Liouville (RL)分數階微分、Caputo分數階微分和Risez分數階微分[16]。其中，Caputo分數階微分在理論求解方面能夠消除常數項對解的影響。因此，下面采用Caputo分數階導數進行具體研究。Caputo分數階導數在區間0

(1)

根據式(1)可得到函數f(t)分別為常函數和冪函數時的通解為：

(a)f(t)=const時：

(2)

(b)f(t)=tx時：

(3)

2 分數階損傷蠕變模型

2.1 損傷描述

損傷力學是通過描述材料的細觀機理來反映其宏觀響應，是一種從細觀唯象理論出發，根據材料自身特性來刻畫結構變形和破壞機理的定量分析方法。周宏偉等[12]提出采用體積應變和波速比來衡量損傷因子的大小。Liao等[8]進一步通過最終應變和初始應變來反映溫凍土的損傷因子。此外，李鑫等[7]通過損傷因子來表征凍土的弱化效應，以此來預測凍土的加速蠕變行為。

由于應力水平和土顆粒與冰晶體的相互作用，使凍土在觸發損傷之后仍然具有繼續承載的能力，并對觸發損傷的應力、時間產生了明顯的影響[7]。因此，對凍土的蠕變曲線(應變ε-時間t曲線)進行分析，具體如圖1所示。其中，Ⅰ代表凍土蠕變最基本的形式，且在此應力作用下不發生黏彈性損傷和蠕變破壞；Ⅱ表明凍土在經歷了衰減蠕變后，發生損傷穩定蠕變和損傷加速蠕變行為，即認為損傷開始于穩定蠕變的初始時刻，并且整個損傷過程歷時較長，因此將其定義為漸進損傷；Ⅲ為應力水平明顯大于損傷閥值，幾乎不經歷穩定蠕變階段就發生蠕變破壞，且整個蠕變過程歷時較短，故稱為快速損傷蠕變行為。

圖1 凍土蠕變的三種基本形式

結合應變對損傷機理的描述和圖1的分析結果，將黏塑性損傷因子考慮為損傷初始應變和最終應變的關系,即：

(4)

式中，D為損傷因子，0≤D≤1；ε為應變；εd為損傷初始應變(見圖1)。根據式(4)和蠕變試驗結果得到損傷因子D與時間te的相互關系。因此，可將式(4)表示為：

D=1-e-αte，

(5)

式中，α為D-te曲線的擬合參數；te為應變為εd時的時間，即te=t-td。

2.2 黏塑性行為

過應力理論是在Malvern[17]一維模型的基礎上發展起來的。在一定程度上，過應力理論也可理解為Bingham模型在三維應力條件下的推廣形式。因此，過應力理論不僅有較為完善的理論基礎，而且在描述凍土蠕變方面有著較為明顯的優勢。Zhou等[13]結合傳統過應力理論和Riesz-Caputo分數階導數獲得了彈黏塑性模型的解析解。因此，在推廣過應力理論的基礎上，考慮分數階導數和損傷效應后，可得到一個黏塑性體的分數階微分關系：

(6)

對式(6)考慮Laplace變換和卷積運算后，得到一般應力狀態下黏塑性體的本構關系：

(7)

式中，p為平均主應力，且p=1/3(σ1+σ2+σ3);Pa為大氣壓力。

此處假設黏塑性變形存在于整個加載過程中，并認為黏塑性損傷只有在應力觸動黏塑性損傷閥值后，才會發生黏塑性破壞。也就是說，當應力小于屈服應力時(f=0)，α=0；當蠕變應力不小于屈服應力時(f>0)，α∈(0,1]。

根據凍土屈服準則和強度效應的研究[18]，表明凍土的臨界狀態線不是直線，而是一條拋物曲線。因此，Liao等[8]根據試驗結果，給出了一個具有一般性的表達式。因此，為了簡化模型，仍然使用拋物線函數來描述凍土的屈服過程：

(8)

式中，ω為拋物線函數與平均應力軸左側交點的橫坐標；φ為內摩擦角；c為黏聚力;q為剪應力;f為屈服函數。

結合屈服方程可得到：

(9)

fq=q+A2，

(10)

根據Caputo分數階導數的運算法則(式(2)和式(3))得到：

(11)

(12)

根據非正交流動法則與式(11)、式(12)可獲得黏塑性流動方向：

(13)

圖2為非正交黏塑性流動法則在子午面上的演化規律。從圖中發現黏塑性流動方向在階次r∈(0,1)時表現出了明顯的非正交性。此外，隨著階次r的減小，非正交性逐漸趨于明顯。值得注意的是，式(13)在沒有引入新勢函數的基礎上就實現了黏塑性流動的非正交特性，從而消除了正交流動在土體的限制[19]。

圖2 非正交黏塑性流動在子午面上的演化

2.3 蠕變模型

為了實現元件模型在描述凍土加載過程中的蠕變變形，采用彈簧體、Abel阻尼器和黏塑性體組建了一種分數階損傷蠕變模型，模型組成具體如圖3所示。

圖3 分數階彈黏塑性蠕變模型示意圖

根據圖3所示，可得到分數階蠕變模型的應力應變關系分別為：

ε=εe+εv+εvp，σ=σe=σv=σvp,

(14)

式中,ε,εe,εv，εvp分別為總應變、彈性應變、黏性應變和黏塑性應變；σ為應力，上標“e”，“v”,“vp”分別表示彈性、黏性和黏塑性性質。

因此，一般應力狀態下凍土的分數階損傷蠕變模型可表示為：

(15)

(16)

式中，G為剪切模量；γa為Abel阻尼器的黏度系數；sij為偏應力張量，sij=σij-pδij；εij(t)為總應變張量。

3 模型驗證

3.1 模型參數分析

模型參數中包含彈性參數、黏性體參數和黏塑性參數。由于彈性參數和黏度系數的普遍性，這里主要討論黏塑性參數對模型計算結果的貢獻。黏塑性體主要包含β，α，te和γ這4個參數。

除此之外，發現采用應變能夠囊括應力、溫度、結構等因素變化對損傷效應的貢獻。所以使用損傷初始應變εd與最終應變的比值來刻畫黏塑性損傷因子，見式(4)。損傷系數α是此類模型反映材料損傷效應的參數之一,是通過擬合損傷因子與時間曲線確定的。由圖4(b)發現黏塑性損傷效應隨著α的增大而快速增大，其中，α=0表現為無損蠕變形式。時間te是除系數α之外的另一個損傷參數,表示損傷時間,te通過時間t和td的差值來計算。其中，td是損傷的觸發時間(見圖1)，損傷觸發時間td對蠕變曲線的影響情況如圖4(c)所示。從圖中發現隨著td的增大，蠕變速率與應變都有少量的降低。故td會影響蠕變的發展速率與應變量級，但不會對其產生本質影響。

圖4 黏塑性體的參數分析

分數階階次r作為分數階模型最主要的一個參數,階次對蠕變曲線的影響情況如圖4(d)所示。從圖中發現在階次r∈(0,1)時可以很好地刻畫材料的非線性蠕變行為。但是，在階次r≥1時應變與時間曲線會出現線性蠕變和蠕變破壞行為。然而，蠕變破壞行為的出現似乎是一個好的結果。

3.2 模型計算

在凍土中，高溫凍土作為一種存在范圍廣、周期長且蠕變效應明顯的特殊凍土，對公路工程和鐵道工程產生了一系列的不良影響。因此，下面分別采用模型模擬高溫凍土的蠕變試驗結果。其中，富含冰凍結沙土在-1 ℃的恒溫下完成加載,溫凍淤泥的試驗溫度為-1.5 ℃。詳細的土性參數和試驗結果分別見Yang等[10]和Liao等[8]。

圖5為模型計算結果與試驗結果的對比情況。其中，圖5(a)為模型對富含冰凍結沙土一維蠕變試驗結果的模擬情況。富含冰凍結沙土的彈性模量和屈服極限分別為E=90.9 MPa和σs=0.568 MPa,其他參數如表1所示。從圖中發現，模型能夠較好地模擬富含冰凍結沙土的穩定蠕變、漸進損傷蠕變和快速損傷蠕變行為。圖5(b)～圖5(d)為模型對溫凍淤泥三軸不排水剪切試驗結果的模擬情況。其中，溫凍淤泥在圍壓為0.5，1.0 MPa和2.0 MPa下的剪切模量G分別為18.13，20.72 MPa和21.13 MPa。可明顯看出，模型計算結果與試驗結果具有較高的一致性。此外，相較于傳統的元件模型而言，新模型是在推廣損傷理論和黏塑性體的基礎上提出的，其表征黏塑性行為也存在于屈服面內，且損傷開始于初始穩定蠕變階段。

圖5 模型對溫凍土蠕變行為的模擬情況

同時，模型也能較好地模擬溫凍淤泥的三維蠕變行為，并且可以刻畫土體在屈服之前的黏塑性行為和屈服之后的損傷蠕變行為。此外，根據模型和試驗結果的對比，還發現無損蠕變過程包含明顯的彈性、黏性變形和少量的黏塑性變形。黏塑性變形會隨著應力水平的增加而增加,具體如表1參數β所示,即β會隨著應力水平的增加而增加[20]。在土體損傷蠕變過程中，黏塑性變形會隨著損傷程度的增加而增加,該現象正好與試驗結果相吻合。

表1 溫凍土的模擬參數

根據模型對凍土蠕變行為的模擬情況，可以肯定模型在模擬凍土蠕變行為的適用性和有效性。具體為，發現推廣后的黏塑性損傷模型可描述凍土屈服前后的黏塑性蠕變行為和損傷效性。除此之外，使用分數階導數表示了凍土黏塑性流動的非正交特性。因此，認為推廣的分數階損傷蠕變模型消除了傳統塑性體純粹塑性的限制。

4 結論

(1)依托試驗結果對凍土損傷蠕變行為進行了分析，將損傷蠕變分為漸進損傷和快速損傷蠕變。在此基礎上，將其進行推廣，并根據試驗結果發現黏塑性損傷模型能夠很好地描述凍土的蠕變破壞行為。

(2)推廣的黏塑性理論考慮了凍土在整個加載過程中的黏塑性行為，消除了塑性體純粹塑性的限制。此外，通過分數階導數，可直接獲得一個非正交黏塑性流動而不需要引入勢函數。

(3)根據試驗結果和模型計算結果的對比情況，發現模型計算結果與試驗結果吻合較好。