慕煥東 鄧亞虹 赫 佳 李 麗

(①長安大學地質(zhì)工程與測繪學院，西安 710054，中國)

(②礦山地質(zhì)災害成災機理與防控重點實驗室，西安 710054，中國)

(③甘肅省交通科學研究院有限公司，蘭州 730050，中國)

0 引 言

黃土廣泛分布于我國西北地區(qū)，其在荷載的長期作用下可能發(fā)生蠕變和漸進破壞等流變現(xiàn)象(如西安地鐵開挖的黃土隧道以及巖土工程中常見的邊坡、擋土墻、地基與路基漸進破壞)，進而影響工程建構(gòu)筑物的長期穩(wěn)定性。因此，分析黃土在長期荷載作用下的應力-應變關系進而建立考慮流變特性的黃土本構(gòu)模型對黃土地區(qū)工程設計具有重要的理論與工程實際意義。

學者們對巖土材料流變特性及本構(gòu)模型的研究多從元件蠕變模型和經(jīng)驗蠕變模型兩方面考慮。經(jīng)驗蠕變模型研究較為廣泛，如陳昌富等(2008，2019)將固結(jié)理論與Kelvin蠕變模型相結(jié)合，建立了考慮固結(jié)-蠕變耦合作用的蠕變模型可合理評價軟土或紅黏土蠕變特性。劉忠玉等(2019)引入考慮時間效應的統(tǒng)一硬化(UH)本構(gòu)模型描述飽和黏土的黏彈塑性并與一維流變固結(jié)試驗結(jié)果對比驗證其有效性。蘇白燕等(2015)對重慶武隆雞尾山滑坡滑帶(炭質(zhì)泥質(zhì)灰?guī)r)進行了巖石飽水直接剪切流變試驗，得出滑帶炭質(zhì)泥質(zhì)灰?guī)r的剪應力-剪切位移曲線及其長期強度。李佐良等(2013)對德州松軟土進行不同圍壓下的三軸固結(jié)不排水蠕變試驗，建立了其分數(shù)階蠕變本構(gòu)模型。

元件模型以Maxwell、Kelvin、Burgers最為經(jīng)典。對于黃土來說，目前大多采用由基本元件模型組成的整數(shù)階模型，林斌等(2010)基于不同含水量Q3黃土的蠕變特征分析，將一個Hooke體與多個Kelvin體串聯(lián)，提出了黃土損傷與流變耦合作用的本構(gòu)模型，該模型可以準確描述黃土的突發(fā)性破壞特征；王鵬程等(2014)對陜西省涇陽縣原狀黃土進行室內(nèi)三軸固結(jié)排水蠕變試驗，以此建立了黏彈塑性九元件蠕變模型，該模型通過1個K體和2個V/K體來描述原狀黃土的線性流變特性，但該模型僅可以較好地描述該地區(qū)黃土的三軸線性蠕變特性，不具備說明問題的普遍性。與此同時，現(xiàn)有研究表明在實際中為了提高模型的預測精度，在整數(shù)階流變本構(gòu)模型中往往通過串聯(lián)多個K體來實現(xiàn)，從而增加了模型的復雜程度及求解參數(shù)的難度。然而，由傳統(tǒng)的線性流變元件串聯(lián)或并聯(lián)組合而成的本構(gòu)模型只能描述黃土的減速與等速蠕變，無法反映加速蠕變階段(孫均，1999)。

上述學者提出的蠕變模型雖具有很多優(yōu)點，但對蠕變試驗曲線中加速蠕變階段不能全面、精準描述，現(xiàn)有的研究還不能形成合理而又統(tǒng)一的認識，是目前研究中所面臨的根本問題。上述問題主要是考慮加速蠕變階段的非線性特性，為此，越來越多的學者們展開對非線性元件本構(gòu)模型的研究。眾所周知，由于含水率的不同，土體的軟硬狀態(tài)不同，因此可以認為土體是一種介于理想固體與理想流體之間的材料。基于這樣的前提，殷德順等(2007)基于分數(shù)階微積分原理，改造出介于理想固體與理想流體之間的軟體元件，建立了改進的分數(shù)階蠕變本構(gòu)模型，與整數(shù)階模型相比，該模型可以更好地描述土體的非線性流變特性。在分數(shù)階流變本構(gòu)模型研究方面，學者們考慮多種類型巖土材料，對其本構(gòu)模型進行拓展得到了分數(shù)階的蠕變本構(gòu)模型，并通過室內(nèi)試驗驗證了其有效性。涉及到的巖土材料主要包括軟黏土、凍土、軟土、泥巖、軟巖、鹽巖等(孫海忠等，2007；何利軍等，2011；宋勇軍等，2013；吳斐等，2014；閆云明等，2017；蔣秀姿等，2018；孫凱等，2018；倪靜等，2019；孫萍萍等，2020)。上述對于分數(shù)階蠕變本構(gòu)模型的研究具有非常重要的理論意義，但目前的研究對黃土材料研究較少，而黃土材料廣泛分布于我國西北地區(qū)，承載的工程建構(gòu)筑物極其廣泛，同時已有研究表明在荷載的長期作用下，黃土所表現(xiàn)出的非線性蠕變特性成為影響其安全的重要因素，尤其是其蠕變變形和長期強度，往往是致使建構(gòu)筑物變形失穩(wěn)亦或強度破壞的主導因素之一。

關于黃土流變特性的研究，本課題組已通過三軸循環(huán)加卸載蠕變試驗，基于統(tǒng)一流變力學模型理論，建立了由西原體與村山體串聯(lián)組合而成的改進西原模型(李麗，2013；鄧亞虹等，2015；李飛霞，2016；李麗等，2017，2018)。該模型能較好地描述Q3黃土的減速蠕變與等速蠕變階段，但對于黃土流變試驗曲線段(減速與加速蠕變段)無法準確描述，尤其是對不同于改進西原模型線性元件的非線性特性的準確描述，其效果有待進一步提高。

鑒于此，本文在課題組前期研究的基礎上，以Q3黃土為研究對象，基于分數(shù)階導數(shù)理論，考慮蠕變特性中等速蠕變及加速蠕變的非線性特征，采用介于理想固體及理想流體之間的Abel黏壺元件代替改進西原模型中的理想固體牛頓黏壺元件，推導建立了分數(shù)階改進西原模型。同時考慮原狀及重塑兩種類型黃土試樣，通過三軸分級循環(huán)加-卸載流變試驗驗證模型的有效性，該研究為實際工程中解決黃土流變問題提供了重要的理論依據(jù)。

1 分數(shù)階導數(shù)基本原理

1.1 基于Riemann-Liouville分數(shù)階微積分理論

相比于整數(shù)階模型，分數(shù)階模型的發(fā)展直至1990年才被廣泛重視，其主要包括了Caputo、Grunwald-Letnikov和Riemann-Liouville(R-L)3種基本形式。上述3種基本形式中：一般通過引入積分定義即可以對微積分的計算過程進行簡化，因此本文采用了R-L積分定義的形式。Riemann-Liouville分數(shù)階微積分的定義為：設f在(0，+∞)上逐段連續(xù)，且在[0，+∞]的任何有限子區(qū)間上可積，對t>0，Re(γ)>0，稱：

(1)

對上述分數(shù)階積分進行逆運算，可以得到Riemann-Liouville分數(shù)階導數(shù)的基本方程為式(2)所示，具體可描述為：設f∈C，ν>0，m是大于γ的最小整數(shù)，記ν=m-γ>0，則可得：

(2)

式(2)稱為函數(shù)f(t)的γ階Riemann-Liouville分數(shù)階導數(shù)(Kilbas et al.，2006)。

1.2 基于分數(shù)階導數(shù)的流變元件

傳統(tǒng)的Newton元件代表著理想流體的含義，Hooke元件代表理想固體，由分數(shù)階導數(shù)建立的流變元件即Abel元件(也被稱為Scott-Blair元件)可以理解為是介于理想固體與理想流體之間的一種本構(gòu)關系。根據(jù)模型中黏滯系數(shù)是否隨時間發(fā)生變化，可將Abel元件劃分為不隨時間變化的常系數(shù)Abel黏壺及隨時間變化的變系數(shù)(主要為考慮損傷特性)Abel黏壺兩種形式(周宏偉等，2012)。常系數(shù)Abel黏壺的本構(gòu)關系式如式(3)所示：

(3)

當σ(t)為常數(shù)時，對式(3)兩側(cè)進行分數(shù)階積分，即可得到常系數(shù)Abel黏壺描述的流變關系式如式(4)所示：

(4)

在流變過程中，黃土的黏性系數(shù)ηγ不再恒為常數(shù)，引入損傷變量D來描述ηγ的劣化，其理論關系式為式(5)：

D=1-e-αt(0≤D<1)

(5)

式中：α為與黃土材料性質(zhì)相關的系數(shù)。

將式(5)代入式(3)可得到考慮損傷變量的變系數(shù)Abel黏壺的本構(gòu)關系為式(6)：

(6)

當σ(t)為常數(shù)時，對式(6)兩側(cè)進行分數(shù)階積分，即可得到變系數(shù)Abel黏壺描述的流變關系如式(7)所示：

(7)

2 流變本構(gòu)模型的建立

2.1 整數(shù)階改進西原模型

黃土的典型流變過程一般分為3個階段，包括減速流變階段、等速流變階段和加速流變階段。夏才初等(2008)根據(jù)不同應力水平下的加卸載蠕變試驗曲線，將彈性、黏性、彈塑性、黏彈性、黏塑性、黏彈塑性力學模型進行不同的串聯(lián)組合，建立了15個流變力學模型，提出了統(tǒng)一流變力學模型理論。鄧亞虹等(2015)基于統(tǒng)一流變力學模型理論，在三軸循環(huán)加卸載流變試驗的基礎上，建立由西原體與村山體串聯(lián)組合而成的整數(shù)階改進西原模型，如圖1所示。

根據(jù)上述所建立的整數(shù)階改進西原模型，考慮西原體中黏塑性體與村山體中黏彈塑性體兩者之間屈服應力的大小關系，可以建立如圖2所示的不同應力水平下的蠕變模型特征曲線。

由圖2所示的蠕變模型特征曲線可知：

(1)在加載瞬間，各應力水平作用下均有瞬時彈性變形，反映出該本構(gòu)模型中應串聯(lián)有單獨的彈簧元件。

(2)在低應力水平作用下，蠕變曲線為減速蠕變類型，卸載后其變形基本能完全回彈，說明具有黏彈性體元件。

(3)在較高應力水平作用下，卸載后蠕變變形有且僅有部分回彈，可以說明此時模型中必然存在黏彈塑性體元件。與此同時，蠕變試驗結(jié)果表明減速蠕變的變形大于卸載后產(chǎn)生的滯后回彈變形，由此可以說明模型中應同時含有黏彈性體元件和黏彈塑性體元件。

(4)在高應力水平作用下，蠕變試驗曲線出現(xiàn)了等速蠕變現(xiàn)象，且最終出現(xiàn)了加速蠕變現(xiàn)象，可見在低應力作用下模型的元件失去了作用，只有當應力達到一定值時才起作用，該元件稱為黏塑性體。

2.2 基于Abel黏壺的分數(shù)階改進西原模型

通過對蠕變模型特征曲線分析發(fā)現(xiàn)，本課題組已經(jīng)建立的如圖1所示的改進西原模型可以較全面地反映Q3黃土在三軸循環(huán)加卸載流變試驗條件下的變形特性，描述其蠕變變形規(guī)律。但深入分析發(fā)現(xiàn)，無論是彈簧原件、黏彈性體元件、黏彈塑性體元件、黏塑性體元件都具有線性元件模型的特征。而黃土是一種典型的非線性土體材料，為了更為有效、真實地模擬黃土的非線性變形特征，在模型中考慮將代表理想流體的牛頓黏壺替換為介于理想固體與理想流體之間的Abel黏壺，建立考慮非線性特性的蠕變本構(gòu)模型。

進一步分析Q3黃土的蠕變曲線，可以總結(jié)出如下特征，在3種應力(低應力、較高應力、高應力)水平作用下呈現(xiàn)出不同的變化特征。具體可描述為低應力水平時的減速蠕變、較高應力水平的減速蠕變、高應力水平的等速和加速蠕變。在第1種情況下，模型中的黏彈性體元件及黏彈塑性體元件的黏滯系數(shù)不隨時間的變化而變化；第2種情況與第1種情況具有相同的變化特征；第3種情況的等速與加速蠕變特性需要黏塑性體來體現(xiàn)，黏滯系數(shù)隨時間的變化而產(chǎn)生變化。因此，通過分析可知，在低應力及較高應力水平作用下，因黏滯系數(shù)為常數(shù)，因此選擇常系數(shù)的Abel黏壺，而在高應力水平作用下，因黏滯系數(shù)不是常數(shù)，因此選擇變系數(shù)的Abel黏壺，由此可構(gòu)建基于Abel黏壺的分數(shù)階改進西原模型如圖3所示。

由圖3可知，分數(shù)階改進西原模型由彈性體、黏彈性體、黏塑性體和黏彈塑性體組成。其本構(gòu)方程按應力大小可描述為不同類型，具體如下：

當σ≤σs1和σs2，即低應力時荷載作用時，首先是模型中彈性體元件的彈簧發(fā)生瞬時彈性變形，其次則為黏彈性體元件發(fā)生減速蠕變變形，在加載完成進行卸載作用，則加載產(chǎn)生的變形完全回彈。根據(jù)式(3)、式(4)關系可得本構(gòu)方程為：

(8)

當σs2<σ<σs1，即較高應力時，彈性體、黏彈性體和黏彈塑性體發(fā)生變形，卸載后變形部分回彈，其本構(gòu)方程為：

(9)

當σ≥σs1和σs2，即高應力時，模型整體均發(fā)生變形，蠕變曲線出現(xiàn)加速蠕變階段，根據(jù)式(6)、式(7)關系可得其本構(gòu)方程為：

(10)

綜上所述，建立分數(shù)階改進西原模型的本構(gòu)方程：

(11)

3 流變試驗與模型驗證

3.1 流變試驗

以西安地區(qū)Q3原狀及重塑黃土為研究對象，采用CSS-2901TS三軸流變試驗儀開展三軸循環(huán)加卸載流變試驗。試驗土樣類型為原狀土樣及重塑土樣兩種，取原狀樣時應避免土樣擾動，取樣完成后應立即裝入圓柱狀鐵皮桶內(nèi)并用膠帶密封；制備重塑土樣時，將取回的一定質(zhì)量的土進行風干、碾碎、過篩、配置含水率，含水率配置時需將水均勻噴灑在土料上表面并攪拌，攪拌后裝袋密封并放置保濕皿至少20h以上。試驗時，控制目標含水率及干密度，采用分層夯實法進行夯實。原狀樣與重塑樣采用同一含水率(20%)，試驗控制條件為固結(jié)不排水(CU試驗)，加載方式為分級循環(huán)加、卸載，加載速率為0.001kPa·s-1，卸載速率為0.005kPa·s-1。固結(jié)時間為24h，固結(jié)后進行加卸載蠕變試驗，衡量變形穩(wěn)定的標準為24h內(nèi)變形量小于0.05mm。試驗完成后將所采集的試驗數(shù)據(jù)按照Boltzmann原理進行處理后，得到了Q3原狀與重塑黃土的蠕變試驗曲線，如圖4所示。

由圖4可知，無論是原狀黃土亦或重塑黃土，隨荷載的作用大小表現(xiàn)出不同的蠕變特性：

(1)在低應力水平階段，蠕變試驗產(chǎn)生的變形較小，如原狀黃土荷載為43.33kPa時，其蠕變變形量為1%，隨著時間的增加蠕變變形量趨于穩(wěn)定，該階段只發(fā)生減速蠕變。

(2)在低應力至較高應力水平階段，蠕變試驗產(chǎn)生的變形有所增加，但增幅較小，相比于低應力階段，其蠕變變形不能完全回彈，說明該階段沒有單獨的黏壺存在。

(3)在高應力水平階段，蠕變試驗產(chǎn)生的變形經(jīng)歷了等速蠕變，并最終出現(xiàn)加速蠕變，可見低應力水平減速蠕變階段的元件模型失去了作用，進而表明模型中存在有黏彈塑性體。

(4)對比原狀黃土和重塑黃土的蠕變試驗曲線，發(fā)現(xiàn)原狀黃土的變形量小于重塑黃土的變形量，變形量差別由1%～7%到10%～40%。

3.2 模型參數(shù)擬合與驗證

為了驗證本文所建立的分數(shù)階本構(gòu)模型的正確性，將其與3.1節(jié)的試驗數(shù)據(jù)繪制在一張圖中，得到不同固結(jié)壓力下原狀黃土和重塑黃土的擬合曲線(圖5、圖6)。從圖可知，模型擬合曲線與試驗數(shù)據(jù)吻合較好，其相關系數(shù)達到了0.99，可見本文所建立的基于Abel黏壺的分數(shù)階改進西原模型本構(gòu)模型可以完整地描述Q3黃土流變的3個階段。

為進一步分析本文所建立的分數(shù)階本構(gòu)模型與鄧亞虹等(2015)所建立的改進西原模型之間所存在差異，將兩者繪制在同一張圖中(圖7)進行分析。由圖7可知，在加載初期(時間較小)的減速蠕變階段，分數(shù)階本構(gòu)模型蠕變應變變化曲線更為圓滑，尤其是重塑黃土體現(xiàn)更為明顯；在應力卸載后的蠕變變形階段，同樣呈現(xiàn)出該現(xiàn)象，說明本文所建立的分數(shù)階本構(gòu)模型可以較好地模擬黃土的流變特性。

在模型驗證的基礎上，需要對模型中各參數(shù)進行求解。最小二乘法作為數(shù)據(jù)處理中確定模型參數(shù)最常用的方法之一，其包含的列文伯格-馬夸爾特法(Levenberg-Marquardt，亦稱之為LM算法)是最小二乘法算法中的一種，其是利用梯度法求最大(小)值來解決曲線最小二乘擬合問題(Madsen et al.，2004；朱杰兵等，2010；何志磊等，2016)，該方法具有梯度法和牛頓法的優(yōu)點，收斂速度更快。

表1 Q3原狀黃土模型參數(shù)擬合結(jié)果Table 1 List of Q3 original loess model parameter fitting results

表2 Q3重塑黃土模型參數(shù)擬合結(jié)果Table 2 List of Q3 remodeling loess model parameter fitting results

4 參數(shù)敏感性分析

根據(jù)表1擬合結(jié)果，將擬合得到第3組參數(shù)代入式(11)，僅改變求導階次γ，可得到不同γ對應的蠕變曲線(圖8)。

由圖8可知，求導階次對蠕變特性具有一定的影響，具體可表現(xiàn)為求導階次為0.1次時，其應變小于2.5%，隨著求導階次的逐漸增加，其應變逐漸增加，當求導階次為0.9時，其應變近2.9%，表明蠕變變形在逐漸增加，但瞬時彈性變形速率則隨著求導階次的增加逐漸減小。

同時，考慮與黏性系數(shù)有關的材料參數(shù)α對蠕變特性的影響，繪制不同材料參數(shù)α下蠕變應變與時間關系曲線如圖9所示。由圖可知，材料參數(shù)為0.1時，其蠕變應變最小，隨著材料參數(shù)的增加，其蠕變應變逐漸增加，但增加的幅值較小，其瞬時彈性變形基本一致，區(qū)別在于最終蠕變應變的大小具有一定差異，材料參數(shù)由0.1變化至0.9時，應變差異近10%。

5 結(jié) 論

(1)基于分數(shù)階導數(shù)理論，考慮蠕變特性中等速蠕變及加速蠕變的非線性特征，引入一種介于理想固體與理想流體之間的分數(shù)階常系數(shù)與變系數(shù)Abel黏壺，建立了分數(shù)階改進西原模型并推導出本構(gòu)方程。

(2)基于Q3原狀與重塑黃土的三軸循環(huán)加卸載流變試驗數(shù)據(jù)，對模型進行了對比分析。結(jié)果表明分數(shù)階改進西原模型既能實現(xiàn)對Q3黃土減速蠕變、等速蠕變和加速蠕變3個流變階段的模擬，可以彌補整數(shù)階改進西原模型無法描述加速蠕變階段，所建立的方程能很好地預測其蠕變變形特性及變形規(guī)律，所得模型參數(shù)可以直接用于Q3原狀與重塑黃土的蠕變數(shù)值分析。

(3)通過對模型參數(shù)進行敏感性分析，發(fā)現(xiàn)求導階次和與黏性系數(shù)有關的材料參數(shù)增加則蠕變變形均逐漸增加，但求導階次增加瞬時彈性變形速率逐漸減小，與黏性系數(shù)有關的材料參數(shù)瞬時彈性變形基本不變。