改進西原蠕變模型研究及二次開發(fā)計算*

王寶迪，宿利平，劉 洋

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

0 引言

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

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

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

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

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

常用的經(jīng)驗模型有冪函數(shù)模型，對數(shù)型模型和指數(shù)型模型等，其基本表達式如下：

1)冪函數(shù)模型

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

ε(t)=a+btc

(1)

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

(2)

式中：a,b,c是經(jīng)驗常數(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)驗常數(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)驗常數(shù)。

1.2 擬合結(jié)果分析

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

(9)

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

表1 軸壓為83.49 kPa時3種經(jīng)驗蠕變模型的比較結(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)驗蠕變模型的比較結(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ù)模型雖然在低應力水平和高應力水平下都有最高的擬合度，但分別對于另1種應力水平的擬合度較差；冪函數(shù)模型在低應力水平的擬合度與指數(shù)模型的擬合度相近，但在高應力水平下，擬合度明顯優(yōu)于指數(shù)函數(shù)的擬合度。

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

冪函數(shù)經(jīng)驗模型對黏土的三軸蠕變曲線在低應力水平下和高應力水平下都有著較高的擬合度。本文在冪函數(shù)經(jīng)驗模型的基礎上，提出冪函數(shù)非線性元件，如圖1所示，其本構如式(10)所示：

(10)

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

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

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

(11)

(12)

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

2 改進西原蠕變模型

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體能夠描述土體在外力作用下達到屈服應力所產(chǎn)生的黏塑性變形。一維西原模型蠕變本構方程如式(13)所示：

(13)

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

(14)

2.2 改進西原蠕變模型

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

圖3 改進西原蠕變模型Fig.3 Modified Nishihara creep model

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

(15)

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

(16)

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

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

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

2.4 改進西原模型適用性驗證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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