基于分?jǐn)?shù)階導(dǎo)數(shù)理論的Q3黃土流變本構(gòu)模型研究*

慕煥東 鄧亞虹 赫 佳 李 麗

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

(②礦山地質(zhì)災(zāi)害成災(zāi)機(jī)理與防控重點(diǎn)實(shí)驗(yàn)室，西安 710054，中國)

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

0 引 言

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

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

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

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

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

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

1 分?jǐn)?shù)階導(dǎo)數(shù)基本原理

1.1 基于Riemann-Liouville分?jǐn)?shù)階微積分理論

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

(1)

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

(2)

式(2)稱為函數(shù)f(t)的γ階Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)(Kilbas et al.，2006)。

1.2 基于分?jǐn)?shù)階導(dǎo)數(shù)的流變元件

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

(3)

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

(4)

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

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

(5)

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

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

(6)

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

(7)

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

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

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

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

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

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

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

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

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

2.2 基于Abel黏壺的分?jǐn)?shù)階改進(jìn)西原模型

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

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

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

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

(8)

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

(9)

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

(10)

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

(11)

3 流變試驗(yàn)與模型驗(yàn)證

3.1 流變試驗(yàn)

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

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

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

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

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

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

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

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

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

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

表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)，僅改變求導(dǎo)階次γ，可得到不同γ對應(yīng)的蠕變曲線(圖8)。

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

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

5 結(jié) 論

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

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

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