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

王寶迪，宿利平，劉 洋

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

0 引言

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

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

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

1 經驗蠕變模型

1.1 3種常用的經驗蠕變模型

常用的經驗模型有冪函數模型，對數型模型和指數型模型等，其基本表達式如下：

1)冪函數模型

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

ε(t)=a+btc

(1)

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

(2)

式中：a,b,c是經驗常數。

2)對數模型

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

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

(3)

ε(t)=bln(1+ct)

(4)

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

(5)

式中：a,b,c,d是經驗常數。

3)指數模型

指數模型的基本形式可以寫成如式(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是經驗常數。

1.2 擬合結果分析

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

(9)

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

表1 軸壓為83.49 kPa時3種經驗蠕變模型的比較結果Table 1 Comparison results of three empirical creep models when the axial pressure is 83.49 kPa

表2 軸壓為333.96 kPa時3種經驗蠕變模型的比較結果Table 2 Comparison results of three empirical creep models when the axial pressure is 333.96 kPa

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

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

1.3 冪函數非線性元件

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

(10)

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

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

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

(11)

(12)

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

2 改進西原蠕變模型

2.1 傳統一維西原蠕變模型

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

圖2 經典西原蠕變模型Fig.2 Classic Nishihara Creep Model

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

(13)

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

(14)

2.2 改進西原蠕變模型

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

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

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

(15)

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

(16)

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

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

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

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

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

表3 經典西原模型擬合所得計算參數Table 3 Calculated parameters obtained from the fitting of the classical Nishihara model

表4 改進西原模型擬合所得計算參數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)所示，經典西原模型和改進的西原模型擬合精度都比較高，決定系數都達到0.98以上，其中改進西原模型的擬合精度相比較經典西原模型的擬合精度更高一些，決定系數達到0.99以上。在軸向應力為250.47 kPa下，如圖5(c)所示，軟土相比較前2個應力水平下的蠕變速率有所增加，經典西原模型的決定系數為0.959 0，改進的西原模型的決定系數為0.982 5，可以看出，改進的西原模型在蠕變速率相對較大的情況下擬合效果更好。在軸向應力為333.96 kPa下，如圖5(d)所示，軟土出現了加速變形階段，經典西原模型的決定系數為0.932 2，改進的西原模型的決定系數為0.991 1，改進的西原模型能夠較好地反映軟土在高應力條件下的瞬時變形、穩態蠕變和加速蠕變過程，整體擬合精度較高。

3 改進西原蠕變模型有限元實現及算例

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

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

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

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

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

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

表6 數值模擬參數Table 6 Numerical simulation parameters

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

4 結論

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

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

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