王寶迪，宿利平，劉 洋

(1.北京科技大學(xué) 土木與資源工程學(xué)院，北京 100083；2.北京市政路橋股份有限公司，北京 100045)

0 引言

流變是巖土體重要的力學(xué)特征之一，包括蠕變、應(yīng)力松弛、長期強(qiáng)度等現(xiàn)象[1-2]。軟黏土具有明顯的時效性,其蠕變特性導(dǎo)致大量工程出現(xiàn)嚴(yán)重的工后沉降問題，如建筑物的下沉及破壞、城市給水及供氣等市政設(shè)施無法正常運(yùn)行、路面高低不平影響交通運(yùn)輸?shù)?。軟土本?gòu)模型的選取是合理預(yù)測軟黏土沉降的關(guān)鍵。目前，關(guān)于巖土體蠕變的研究已取得一定成果，其中，經(jīng)驗(yàn)?zāi)Ｐ秃驮Ｐ蜑槟壳斑\(yùn)用最廣泛的蠕變本構(gòu)模型[3]。經(jīng)驗(yàn)?zāi)Ｐ途哂斜磉_(dá)式簡單、模型參數(shù)易獲取和運(yùn)用方便等特點(diǎn)。目前已有不少國內(nèi)外學(xué)者基于大量試驗(yàn)，應(yīng)用經(jīng)驗(yàn)?zāi)Ｐ脱芯坎煌瑧?yīng)力水平下軟土蠕變特性[4-7]。但其研究成果主要基于三軸蠕變試驗(yàn)對不同應(yīng)力水平下蠕變曲線進(jìn)行研究，選用某種特定的函數(shù)形式與蠕變曲線對比分析，而未對不同函數(shù)形式的經(jīng)驗(yàn)蠕變模型進(jìn)行分析評價。

常用的元件模型有Kelvin模型、Merchant模型及西原模型等。其中傳統(tǒng)西原模型在三參量模型H-K的基礎(chǔ)上串聯(lián)了黏塑性元件，能較好地模擬巖土的流變特征。一些學(xué)者為了更準(zhǔn)確地描述巖土體流變規(guī)律，在傳統(tǒng)西原模型上作出改進(jìn)[8-11]。以上研究成果主要對經(jīng)典元件模型做修正，從而改進(jìn)模型在描述巖土體流變階段中的不足，研究對象多為巖體，而對于軟土的研究相對較少。

本文基于經(jīng)驗(yàn)?zāi)Ｐ偷难芯刻岢鰞绾瘮?shù)非線性元件，并在經(jīng)典西原模型基礎(chǔ)上進(jìn)行改進(jìn)，提出改進(jìn)西原蠕變模型。根據(jù)試驗(yàn)結(jié)果[12-13]研究經(jīng)典西原模型與改進(jìn)西原模型描述不同應(yīng)力狀態(tài)下的準(zhǔn)確度及參數(shù)敏感性，基于有限元軟件ABAQUS平臺對改進(jìn)西原蠕變模型進(jìn)行二次開發(fā)，建立數(shù)值模型，將程序模擬值與試驗(yàn)值進(jìn)行對比，驗(yàn)證本文建立的改進(jìn)西原模型及子程序的適用性及有效性。改進(jìn)西原模型及相關(guān)二次開發(fā)可對軟黏土路基進(jìn)行快速、準(zhǔn)確的沉降預(yù)測，以保證對軟土工程設(shè)計、施工、加固、災(zāi)害預(yù)防和防治全工程周期提出合理化建議以達(dá)到工程安全生產(chǎn)及運(yùn)行的目的，為軟土的沉降計算提供新的途徑和方法。

1 經(jīng)驗(yàn)蠕變模型

1.1 3種常用的經(jīng)驗(yàn)蠕變模型

常用的經(jīng)驗(yàn)?zāi)Ｐ陀袃绾瘮?shù)模型，對數(shù)型模型和指數(shù)型模型等，其基本表達(dá)式如下：

1)冪函數(shù)模型

冪函數(shù)模型的基本形式可以寫成如式(1)～(2)所示：

ε(t)=a+btc

(1)

ε(t)=a[(1+bt)c+1]

(2)

式中：a,b,c是經(jīng)驗(yàn)常數(shù)。

2)對數(shù)模型

對數(shù)模型的基本形式可以寫成如式(3)～(5)所示：

ε(t)=a+bln(t+c)+dt

(3)

ε(t)=bln(1+ct)

(4)

ε(t)=a+bln(t+c)

(5)

式中：a,b,c,d是經(jīng)驗(yàn)常數(shù)。

3)指數(shù)模型

指數(shù)模型的基本形式可以寫成如式(6)～(8)所示：

ε(t)=A(1-be-ct)

(6)

ε(t)=a+A(1-e-ct)

(7)

ε(t)=a+A(1-be-ct)

(8)

式中：a,b,c,A是經(jīng)驗(yàn)常數(shù)。

1.2 擬合結(jié)果分析

根據(jù)以上模型，選取成都黏土蠕變試驗(yàn)結(jié)果[12]作為參數(shù)擬合的參照樣本，使用最小二乘法對上述8個模型進(jìn)行參數(shù)擬合，本文以基準(zhǔn)指數(shù)決定系數(shù)(R2)用于評估上述8個模型的性能，如式(9)所示：

(9)

從表1～2可知，對于軸壓為83.49 kPa的蠕變曲線，對數(shù)模型的式(4)的決定系數(shù)最高，擬合度最好，但對于高軸壓333.96 kPa的決定系數(shù)僅為0.039 5，不適合應(yīng)用于高應(yīng)力水平的蠕變曲線計算；對于軸壓為333.96 kPa的蠕變曲線，對數(shù)模型的式(5)的決定系數(shù)最高，擬合度最好，但對于低軸壓83.49 kPa的蠕變曲線無法得到曲線擬合結(jié)果，不適合應(yīng)用于低應(yīng)力水平的蠕變曲線計算。

表1 軸壓為83.49 kPa時3種經(jīng)驗(yàn)蠕變模型的比較結(jié)果Table 1 Comparison results of three empirical creep models when the axial pressure is 83.49 kPa

表2 軸壓為333.96 kPa時3種經(jīng)驗(yàn)蠕變模型的比較結(jié)果Table 2 Comparison results of three empirical creep models when the axial pressure is 333.96 kPa

冪函數(shù)模型和指數(shù)模型在軸壓為83.49 kPa下的擬合度差別較小，且都能大于0.90以上；而對于軸壓為333.96 kPa蠕變曲線，冪函數(shù)模型的擬合度大于0.90，明顯優(yōu)于指數(shù)模型的擬合度。

綜上所述，對數(shù)模型雖然在低應(yīng)力水平和高應(yīng)力水平下都有最高的擬合度，但分別對于另1種應(yīng)力水平的擬合度較差；冪函數(shù)模型在低應(yīng)力水平的擬合度與指數(shù)模型的擬合度相近，但在高應(yīng)力水平下，擬合度明顯優(yōu)于指數(shù)函數(shù)的擬合度。

1.3 冪函數(shù)非線性元件

冪函數(shù)經(jīng)驗(yàn)?zāi)Ｐ蛯︷ね恋娜S蠕變曲線在低應(yīng)力水平下和高應(yīng)力水平下都有著較高的擬合度。本文在冪函數(shù)經(jīng)驗(yàn)?zāi)Ｐ偷幕A(chǔ)上，提出冪函數(shù)非線性元件，如圖1所示，其本構(gòu)如式(10)所示：

(10)

式中：Em′為冪函數(shù)元件簡化前的初始彈性模量；a′，b′，c′為冪函數(shù)元件簡化前的計算參數(shù)。

圖1 冪函數(shù)非線性元件Fig.1 Power function nonlinear element

這里假設(shè)a′不為0，故冪函數(shù)非線性元件本構(gòu)可化簡為如式(11)～(12)所示：

(11)

(12)

式中：Em為冪函數(shù)非線性彈性模量；E0為初始彈性模量；b，c為冪函數(shù)非線性元件的計算參數(shù)。

2 改進(jìn)西原蠕變模型

2.1 傳統(tǒng)一維西原蠕變模型

西原模型由1個Hooke體、1個Kelvin體和1個非線性Bingham體串聯(lián)而成。如圖2所示。

圖2 經(jīng)典西原蠕變模型Fig.2 Classic Nishihara Creep Model

Hooke體能夠描述土體在外力作用下的瞬時變形；Kelvin體將1個Hooke體和1個黏壺元件并聯(lián)，能夠描述土體在外力作用下的黏彈性變形；非線性Bingham體能夠描述土體在外力作用下達(dá)到屈服應(yīng)力所產(chǎn)生的黏塑性變形。一維西原模型蠕變本構(gòu)方程如式(13)所示：

(13)

式中：E0為Hooke體的彈性模量；E1，η1分別為Kelvin體的彈性模量和黏滯系數(shù)；σs，η2分別為Bingham體的屈服應(yīng)力和黏滯系數(shù)。定義如式(14)所示：

(14)

2.2 改進(jìn)西原蠕變模型

傳統(tǒng)的西原模型被廣泛應(yīng)用于描述巖土材料的蠕變特性，能夠描述衰減和穩(wěn)態(tài)蠕變過程，但不能較準(zhǔn)確地反映加速蠕變并且對非線性性質(zhì)的描述還存在不足。為了較清楚地描述不同應(yīng)力水平下黏土蠕變?nèi)^程，本文在原始西原模型的基礎(chǔ)上應(yīng)用上文提出的冪函數(shù)非線性元件替換西原模型原有的彈性元件，并引入非線性黏塑性體(NVPB)模型[14]構(gòu)成改進(jìn)的西原模型，其模型結(jié)構(gòu)示意圖如圖3所示。

圖3 改進(jìn)西原蠕變模型Fig.3 Modified Nishihara creep model

NVPB模型由1個塑性開關(guān)和非線性黏性元件并聯(lián)組成，在恒力σ0下的蠕變方程如式(15)所示：

(15)

式中：ε為在軟土應(yīng)力σ下的應(yīng)變量；ηm為NVPB體黏性系數(shù)；t為蠕變總時間；t0為參考時間。改進(jìn)的一維西原模型蠕變本構(gòu)方程如式(16)所示：

(16)

2.3 改進(jìn)西原模型蠕變參數(shù)敏感性分析

在應(yīng)力狀態(tài)相同的條件下，對式(16)的待定系數(shù)b，c，n進(jìn)行敏感性分析。圖4(a)～圖4(c)分別為系數(shù)b，c，n對改進(jìn)西原模型蠕變曲線的影響。圖4(a)可以看出b的變化對土體初始瞬時應(yīng)變大小影響較為顯著，而對蠕變曲線的最終總應(yīng)變量影響較小，b的值越大蠕變曲線的初始應(yīng)變值越大。圖4(b)可以看出c的變化主要影響蠕變曲線的最終總應(yīng)變量，當(dāng)c值減小并小于零時，蠕變曲線應(yīng)變值明顯增大，當(dāng)c值增大并大于零時，蠕變曲線可以描述土體的負(fù)蠕變現(xiàn)象。圖4(c)可以看出n的取值對蠕變曲線所反映的蠕變變形有著顯著影響：當(dāng)n≤1時，蠕變曲線可描述黏土瞬時彈性變形和衰減蠕變變形；當(dāng)n>1時，蠕變速率逐漸增加，且n越大，蠕變速率越大。

圖4 蠕變參數(shù)對改進(jìn)西原模型蠕變曲線的影響Fig.4 The influence of creep parameters on the creep curve of the improved Nishihara model

2.4 改進(jìn)西原模型適用性驗(yàn)證

為驗(yàn)證本模型的精度和適用性，選取成都黏土[12]蠕變試驗(yàn)結(jié)果作為參數(shù)擬合的參照樣本，該文獻(xiàn)以成都黏土作為研究對象進(jìn)行不同應(yīng)力下的三軸蠕變試驗(yàn)，軟土在低應(yīng)力條件下(軸向應(yīng)力為83.49，166.98，250.47 kPa)其中變形經(jīng)過瞬時蠕變和衰減蠕變后趨于穩(wěn)定，而在高應(yīng)力條件下(軸向應(yīng)力為333.96 kPa)，蠕變變形最終進(jìn)入加速蠕變而使土體發(fā)生破壞。本文采用成都黏土試驗(yàn)蠕變結(jié)果對經(jīng)典西原模型及改進(jìn)的西原模型進(jìn)行擬合分析，擬合得到的計算參數(shù)見表3～4，擬合曲線見圖5。

表3 經(jīng)典西原模型擬合所得計算參數(shù)Table 3 Calculated parameters obtained from the fitting of the classical Nishihara model

表4 改進(jìn)西原模型擬合所得計算參數(shù)Table 4 Improved calculation parameters obtained from Nishihara model fitting

圖5 文獻(xiàn)[12]蠕變試驗(yàn)的擬合曲線Fig.5 Fitting curves of creep tests in Reference [12]

軟土在軸向應(yīng)力為83.49，166.98 kPa下，如圖5(a)～5(b)所示，經(jīng)典西原模型和改進(jìn)的西原模型擬合精度都比較高，決定系數(shù)都達(dá)到0.98以上，其中改進(jìn)西原模型的擬合精度相比較經(jīng)典西原模型的擬合精度更高一些，決定系數(shù)達(dá)到0.99以上。在軸向應(yīng)力為250.47 kPa下，如圖5(c)所示，軟土相比較前2個應(yīng)力水平下的蠕變速率有所增加，經(jīng)典西原模型的決定系數(shù)為0.959 0，改進(jìn)的西原模型的決定系數(shù)為0.982 5，可以看出，改進(jìn)的西原模型在蠕變速率相對較大的情況下擬合效果更好。在軸向應(yīng)力為333.96 kPa下，如圖5(d)所示，軟土出現(xiàn)了加速變形階段，經(jīng)典西原模型的決定系數(shù)為0.932 2，改進(jìn)的西原模型的決定系數(shù)為0.991 1，改進(jìn)的西原模型能夠較好地反映軟土在高應(yīng)力條件下的瞬時變形、穩(wěn)態(tài)蠕變和加速蠕變過程，整體擬合精度較高。

3 改進(jìn)西原蠕變模型有限元實(shí)現(xiàn)及算例

本文利用ABAQUS軟件提供的材料自定義本構(gòu)開發(fā)接口UMAT，對改進(jìn)西原模型進(jìn)行二次開發(fā)。

針對南沙軟土一維蠕變試驗(yàn)結(jié)果[13]運(yùn)用MATLAB軟件的最小二乘法進(jìn)行參數(shù)擬合，參數(shù)擬合結(jié)果及精度如表5所示。由表5可知，得出的模型參數(shù)E0變化范圍不大，其余參數(shù)變化規(guī)律明顯。如圖6 所示，其中E0的變化范圍在870～930 kPa之間，E1，η1，c隨著軸壓的增加而成指數(shù)增大，而b隨著軸壓的增大而減小。

表5 改進(jìn)西原模型擬合文獻(xiàn)[13]所得計算參數(shù)Table 5 Improve the calculation parameters of Nishihara model fitting literature [13]

圖6 蠕變參數(shù)與應(yīng)力的關(guān)系Fig.6 The relationship between creep parameters and stress

為了驗(yàn)證所建模型及二次開發(fā)過程的可行性和合理性，建立與文獻(xiàn)[14]室內(nèi)試驗(yàn)規(guī)格相同的內(nèi)徑為6.18 cm，高度為8 cm的南沙軟土圓柱體數(shù)值模型，荷載選用與原文室內(nèi)試驗(yàn)相同的荷載，模擬試驗(yàn)共為6組，軸向應(yīng)力分別為25，50，100，200，300，400 kPa。

在數(shù)值模擬過程中設(shè)置邊界條件時，對模型底部進(jìn)行完全固定約束，圓柱面進(jìn)行U1和U2方向的位移約束。模型材料參數(shù)選定時，應(yīng)利用多級應(yīng)力水平下試驗(yàn)結(jié)果求出各模型參數(shù)的數(shù)學(xué)期望，并采用一些統(tǒng)計方法和優(yōu)化方法，將會得到模型參數(shù)更為優(yōu)化的值[15]。E0采用擬合值的平均值，其余參數(shù)根據(jù)表6中改進(jìn)西原模型參數(shù)與應(yīng)力關(guān)系選定。在上述步驟均按順序完成后，在提交作業(yè)時，選中編輯好的用戶子程序文件，得到的結(jié)果與試驗(yàn)結(jié)果作比較，如圖7所示?？梢钥闯鲈谳S向應(yīng)力為25，50，100，200，300，400 kPa下，ABAQUS中應(yīng)用改進(jìn)西原蠕變模型UMAT子程序的模擬試驗(yàn)結(jié)果與室內(nèi)試驗(yàn)結(jié)果吻合良好，證明該程序能夠很好地反映一維下軟土瞬時彈性變形和衰減蠕變變形。需要指出的是，本文未研究三維狀態(tài)下的改進(jìn)西原模型參數(shù)與應(yīng)力之間的關(guān)系及比較軟土在三維應(yīng)力狀態(tài)下室內(nèi)試驗(yàn)結(jié)果與模擬試驗(yàn)結(jié)果，這將是進(jìn)一步要研究的內(nèi)容。

表6 數(shù)值模擬參數(shù)Table 6 Numerical simulation parameters

圖7 數(shù)值模擬結(jié)果與蠕變試驗(yàn)結(jié)果對比曲線Fig.7 Comparison curve between numerical simulation results and creep test results

4 結(jié)論

1)分析3種常用的經(jīng)驗(yàn)蠕變模型對軟土蠕變研究的適用性，對比分析發(fā)現(xiàn)冪函數(shù)在軟土處于低應(yīng)力狀態(tài)下產(chǎn)生的瞬時彈性變形和衰減蠕變變形以及處于高應(yīng)力狀態(tài)下產(chǎn)生的加速蠕變變形都有著良好的擬合度，提出冪函數(shù)非線性元件。

2)改進(jìn)后的西原模型能夠較好地描述軟土的非線性蠕變特性及軟土在高應(yīng)力水平下的加速蠕變變形階段。

3)利用南沙軟土的單軸試驗(yàn)結(jié)果與利用二次開發(fā)子程序模型模擬結(jié)果進(jìn)行對比驗(yàn)證，結(jié)果表明有限元計算模型可以反映不同軸壓下的蠕變關(guān)系，研究結(jié)果可進(jìn)一步應(yīng)用于相關(guān)巖土工程蠕變數(shù)值分析，為巖土工程蠕變分析多樣化提供1種可行的方案。