路基壓實粉質黏土動態回彈模量的四參數預估模型及有限元實現

董 城, 劉文劼, 李 亮

(1.湖南省交通科學研究院有限公司,長沙 410015;2.中南大學 土木工程學院,長沙 410075)

長期以來，我國在路基路面結構分析中采用是基于靜力學的彈性層狀體系理論，假定各層結構為勻質材料，每層材料只由回彈模量和泊松比描述[1-2]。該理論在低速、低軸載條件下可以被接受，但為了適應高速、重載公路的迅猛發展，開展基于動力學理論的路基路面結構設計方法成為了必然。在動力學框架下，結構內部的應力狀態時刻發生變化，因此，采用受應力狀態影響的動態回彈模量作為路基力學特性表征參數將能更好地與動力學分析相結合，更準確地反映路基實際變形情況[3]。

國內外許多學者和研究機構在深入分析應力狀態對動態回彈模量影響規律的基礎上，陸續建立了一些動態回彈模量預估模型以嘗試定量描述路基土的變形行為，其中，考慮剪切影響的模型包括雙線性模型、冪指數模型，考慮側限影響的模型包括圍壓模型、K-θ模型，但這兩類模型的主要缺點在于僅考慮了體應力或偏應力對路基土動態回彈模量的影響，而實際上大部分土與粒料的動態回彈模量是體應力、偏應力兩者的函數，因此綜合考慮體應力和偏應力的影響的復合類模型可以更真實、全面地反映材料的力學性狀，Uzan等[4-5]都相繼提出了復合類預估模型，美國于2004年在NCHRP 1-37A項目中提出公路路基設計的動態回彈模量模型(簡稱N37A模型)[6]，消除了模量不定值和量綱不一致問題，目前被工程界廣泛接受。后續許多回彈模量預估模型均是在N37A回彈模型的基礎上進行若干修正[7-8]。

針對具體路基土，建立合適的動態回彈模量預估模型并將其移植入有限元軟件中，是利用動態回彈模量對路基路面進行結構有限元動力學分析的前提，但由于我國路面結構、所處的自然環境以及路基土材料屬性的差異，一些國外成熟的路基土回彈模量預估模型不能直接運用于我國路基路面結構設計，仍需結合我國的地質條件和實驗數據進行改進。本研究以湖南湘東地區路基粉質黏土為例，通過室內動三軸試驗測定不同應力狀態、壓實度、含水率下粉質黏土的動態回彈模量，分析影響粉質黏土動態回彈模量的因素及相關規律，并在N37A模型的基礎上提出含四參數的N37AP1模型，通過多元回歸分析，得到了不同壓實度下粉質黏土對應的N37AP1模型參數。在此基礎上，嚴格推導出N37AP1模型的一致切線剛度矩陣，編寫用戶自定義材料子程序(UMAT)將其移植到有限元軟件ABAQUS 中，通過單個單元的軸壓、圍壓加載以及典型路基路面結構分析，驗證所推導一致性切線剛度矩陣及UMAT編寫的正確性，為粉質黏土路基—路面結構設計提供更真實有效的數值模擬方法，也可供其它種類路基填土的類似研究借鑒參考。

1 動態回彈模量試驗及數據分析

1.1 試驗概況

本文選取湖南湘東地區具有代表性的粉質黏土路基填料作為研究對象，其基本物理性質指標列于表1。

表1 土樣基本物理性質指標

試驗參照我國《公路路基設計規范》(JTG D30—2015)[9]附錄A進行，將粉質黏土在不同含水率(wopt+3%、wopt、wopt-3%)和不同壓實度K(93%、96%)條件下制備試樣。路基土應力加載序列采用國內已有研究成果，動偏應力范圍為30～105 kPa，圍壓范圍為15～60 kPa，動態回彈模量測試中，某一加載級位下荷載循環次數為100 次，記錄各級重復荷載作用下最后5 次循環的回彈變形平均值作為動態回彈模量的計算依據。

1.2 應力狀態與回彈模量的關系

圖1為w=wopt狀態下粉質黏土動態回彈模量隨應力狀態的變化曲線。根據圖1(a)和圖1(b)，粉質黏土動態回彈模量隨著動偏應力的升高而減小，且減小的速率逐漸增大；而根據圖1(c)和圖1(d)，動態回彈模量將隨著體應力的升高而增大，基本呈線性變化或速率逐漸放緩。總體上看，壓實度對動態回彈模量的影響也比較明顯，如當體應力θ=210 kPa且σd= 30 kPa時，相對93%壓實度，壓實度為96%時動態回彈模量提高了16.65%。

(a)

(b)

(c)

(d)

2 動態回彈模量模型改進及評價

2.1 已有模型的改進及回歸分析

在建立動態回彈模量的預估模型時，首先應盡可能準確地按本構定律建立起應力-應變的關系模型，而其他因素的影響則可通過模型參數予以反映。已有研究和本文試驗均表明，路基土的動態回彈模同時受體應力和偏應力的雙重影響，考慮該現象的復合式模型中具有代表性的N37A模型形式為

(1)

對于不同路基土，剪應力對路基土回彈模量的影響程度也是不同的，據此考慮在N37A模型形式基礎上，引入k4項用于調節偏應力對回彈模量的貢獻比例，得到的改進模型(簡稱N37AP1)具體為

(2)

當k4取很小值時，該模型逼近于K-θ模型[10]；當k4取1則成為N37A模型。

分別采用N37AP1模型、N37A模型和UZANP1[11]模型對對本次試驗結果進行了回歸分析，得到了粉質黏土在不同條件下的預估模型參數和對應的決定系數R[12-13]，回歸結果見表2。

表2 預估模型參數回歸結果(k4=0.1)

由表2可見，在不同的含水率和壓實度下，UZANP1模型的決定系數最低，N37A模型次之，而N37AP1模型的決定系數在同條件下都大于其它兩者，對于試驗結果的回歸效果明顯更好。

2.2 不同模型預估結果對比

為了驗證N37AP1模型的預估效果，測試了更多應力狀態下粉質黏土的動態回彈模量，表3給出了w=wopt，K=96%條件下粉質黏土動態回彈模量實測數據和三種模型預估值，總體上改進模型(k4=0.1)的每項誤差值都最小，而且在剪應力較小時誤差都在5%以下，當剪應力增大時誤差略有增加，最大為8.3%。

綜上分析可知，N37AP1模型對最佳含水率下的粉質黏土的動態回彈模量的預估能力明顯更加突出，誤差大體上是N37A模型的50%～60%，是UZANP1模型的20%～45%。N37AP1模型通過調整k4取值的大小，改變偏應力對粉質黏土動態回彈模量的貢獻，進而改變了動態回彈模量的預估結果，使得預估模型能更好的擬合試驗數據。

3 N37AP1一致性切線剛度矩陣的推導

建立動態回彈模量預估模型的主要目的之一為了將其有效應用于有限元計算，本文利用ABAQUS的UMAT接口將N37AP1模型移植入該軟件中。ABAQUS /standard分析模塊的每一步非線性求解都是通過采用切線剛度矩陣進行的，因此要實現N37AP1模型的二次開發就必須求得其全應力狀態下的一致切線剛度矩陣，求解過程如下：

線彈性本構關系可表示為

(3)

式中：E為應變張量；S為應力張量；E為線彈性模量；ε=tr(E)為體應變；v為泊松比；α=ν/(1-2ν)。

類似于線彈性本構關系，非線性彈性本構關系可以利用剛度矩陣Mr寫成

表3 動態回彈模量預測與試驗結果對比(k4=0.1)

(4)

令

(5)

式中：k=k1/(1+ν)；k1,k2,k3,k4為N37AP1模型參數；θ為體應力；Pa為標準大氣壓；τoct為八面體剪應力。

結合式(4)和式(5)，動態回彈模量本構關系可以簡寫為

S=C(θ,τoct)(αεI+E)

(6)

由式(6)可得體應力表達式

(7)

體應力和偏應力可以簡寫為

(8)

(9)

將C(θ,τoct)的表達式代入式(9)，則式(9)可以寫成

(10)

顯然θ,τoct與ζ,γ之間存在函數關系

(11)

因此，C=C(θ,τoct)可以用應變表示為：C=C(ζ,γ)，進而，式(6)可以表示成

S=C(ζ,γ)(αεI+E)

(12)

對式(12)求導可以得到材料的切線剛度矩陣D=?S/?E

(1+αI?I)+(αεI+E)?▽EC

(13)

由式(12)和式(13)可得

(14)

由定義υ=|ε|，可得

(15)

(16)

(17)

由式(5)可得

(18)

由式(10)可得

(19)

(20)

結合式(9)將切線剛度矩陣寫成偏應力的形式

(21)

為方便編程，將切線剛度矩陣寫成如下形式

(22)

按照以上推導結果，在UMAT中以當前應力張量計算出一致切線剛度矩陣(雅可比矩陣)[14]，然后結合當前應變增量得到應力增量，最后根據當前應力和應力增量更新當前的應力，將更新后的應力傳入下一步計算，即可實現預估模型在ABAQUS軟件中的二次開發(見圖2)。

圖2 UMAT實現流程圖

4 預估模型有限元實現驗證

4.1 預估模型的單元驗證

為了驗證所推導一致性切線剛度矩陣及UMAT編寫的正確性，在ABAQUS中先采用單個單元進行驗證，驗證過程中k4取不同值。采用材料參數見表4。

表4 N37AP1模型的材料參數

有限元模型采用單個單元C3D8，約束相鄰三個面的法向位移，在其余三個自由面上施加法向荷載。軸力和側壓力均為漸變形式加載，初始荷載：σ1=0 kPa，σ2=σ3=0 kPa；最終荷載：軸力σ1=100 kPa，σ2=σ3=10 kPa,如圖3所示。

圖3 單元加載示意圖

應變的解析解可表達為

(23)

提取單個單元回彈模量隨軸力變化曲線以及軸向應力應變和側向應力應變曲線，如圖4～圖6所示。

圖4 Mr與σ1的關系曲線

圖5 ε1與σ1關系曲線

圖6 ε3與σ3關系曲線

由圖4～圖6可知，在漸變加載工況下，三組參數所得的應力應變曲線和理論值十分接近，因而推導的一致切向剛度矩陣和編寫的UMAT是正確的。同時，也可以看出k4在改進模型中所起的作用，能調整八面體剪切應力在預估模型中的比重，在相同應力狀態下，k4值越大回彈模量預估值越小，預估應變越大，路基土剪切軟化效果越明顯；反之則回彈模量預估值越大，預估應變越小。

4.2 預估模型的結構驗證

為進一步驗證預估模型，將其移植入了路基—路面有限元模型中，該模型尺寸參照高速公路粉質黏土路基—瀝青路面結構，取 1/4模型進行研究，如圖7所示，模型參數參考湖南省某高速公路的設計參數，如表5所示，其中路基土的材料參數采用w=wopt時粉質黏土的動態回彈模量改進模型參數。車輛荷載形式為單軸雙輪組，軸載大小為100 kN，為了便于設置輪胎接地區域，將該區域等效為與網格形狀一致的矩形[15]，考慮行車道和超車道同時存在車輛的情況。

表5 路基路面結構參數

圖8為在自重作用下路基路面結構路基土回彈模量的分布云圖，由于路基土采用了應力相關的動態回彈模量，路基內不同位置的回彈模量值各不相同。

圖7 整體網格

圖8 自重作用下路基土回彈模量云圖

圖9為路基各點回彈模量在豎向和橫向上的分布曲線。從圖9(a)可以看出，在自重作用下路基土回彈模量隨深度增加而有所增加，最終趨于穩定，路床的回彈模量值在70～100 MPa，路基本體的回彈模量值在60～85 MPa，由于K=0.96和K=0.93路基土的模型參數不同，回彈模量在路床底有了突變，體現了上下層壓實度差異造成的影響。從圖9(b)可以看出，路基的回彈模量由路基中央向兩側衰減，但深度越大，衰減效果越不明顯。

圖10為施加車輛荷載后，由于應力狀態變化形成的回彈模量相對變化率云圖。對于行車道和超車道下方的路基，由于偏應力的明顯增加，回彈模量出現了衰減，最大衰減幅為-1.3%；路床層的車輪間隙區域出現回彈模量的增加，最大增幅約2.0%。由此可以看出，采用動態回彈模量預估模型有效體現了路基各點回彈模量因車輛荷載作用產生的時間變化。

(a)豎向分布

(b)橫向分布

圖10 施加車輛荷載后路基土回彈模量的相對變化率云圖

5 結 論

(1)結合粉質黏土動態回彈模量試驗數據，對已有的動態回彈模量預估模型進行了改進，提出綜合N37A模型、UZANP1模型和K-θ模型的四參數改進模型N37AP1，確定了N37AP1模型對不同含水率和壓實度的粉質黏土的預估參數，通過與實測值的比較，驗證了改進模型對粉質黏土動態回彈模量的預估效果更為理想。

(2)基于廣義虎克定律推導了N37AP1模型的精確一致切線剛度矩陣，通過用戶子程序UMAT成功將該模型嵌入到有限元軟件ABAQUS。

(3)單個單元測試結果表明，k4在改進模型中能調整八面體剪切應力在預估模型中的比重，k4=0.1、k4=1、k4=10時所得的應力應變曲線都和理論值十分接近，在相同應力狀態下，k4值越大回彈模量預估值越小，預估應變越大，路基土剪切軟化效果越明顯。

(4)對典型粉質黏土路基—瀝青混凝土路面結構的分析結果表明，對路基土采用N37AP1模型可以有效反映路基各點在自重作用下回彈模量的空間差異，且可以體現壓實度的影響；同時能反映車輛荷載作用下回彈模量的演變，實現了模量與應力狀態的動態耦合。四參數預估模型在有限元軟件中的嵌入及運用為路基路面結構設計提供了更為真實有效的數值模擬方法。

參 考 文 獻

[1] 羅志剛. 路基與粒料層動態模量參數研究[D].上海:同濟大學, 2007.

[2] 鄧學鈞.路基路面工程[M].北京:人民交通出版社,2010.

[3] 曹海瑩,朱毅,劉云飛,等. 車輛荷載作用下雙層路基層間動應力響應試驗研究[J]. 振動與沖擊, 2017, 36(5): 6-10.

CAO Haiying,ZHU Yi,LIU Yunfei,et al. Testing research on dynamic stress response at interlayer of two-layer roadbed under vehicle load[J]. Journal of Vibration and Shock, 2017, 36(5): 6-10.

[4] LEKARP F, ISACSSON U, DAWSON A. State of the art I: resilient response of unbound aggregates[J]. Journal of Transportation Engineering,2000,126(1):66-75.

[5] LI D, SELIG E T. Resilient modulus for fine-grained subgrade soil[J]. Journal of Geotechnical Engieering, ASCE, 1994,120(6):939-957.

[6] 董城,冷伍明,李志勇.粉土動態回彈模量試驗研究[J].中南大學學報(自然科學版),2012,43(12):4834-4839.

DONG Cheng, LENG Wuming, LI Zhiyong. Dynamic resilient modulus of silt[J]. Journal of Central South University(Science and Technology), 2012,43(12):4834-4839.

[7] 陳聲凱,凌建明,羅志剛. 路基土回彈模量應力依賴性分析及預估模型[J].土木工程學報, 2007, 40(6): 95-99.

CHEN Shengkai, LING Jianming, LUO Zhigang.Stress-dependent characteristics and prediction model of the resilient modulus of subgrade soils[J]. China Civil Engineering Journal, 2007, 40(6): 95-99.

[8] 凌建明,蘇華才,謝華昌,等. 路基土動態回彈模量的試驗研究[J]. 地下空間與工程學報,2010,6(5): 919-925.

LING Jianming, SU Huacai, XIE Huachang, et al. Library research on dynamic resilient modulus of subgrade soil[J]. Chinese Journal of Underground Space and Engineering, 2010,6(5): 919-925.

[9] 中華人民共和國交通運輸部公路路基設計規范:JTG D 30—2015[S].北京:人民交通出版社,2015.

[10] WITCZAK M W,UZAN J. The universal airport pavement design system[R]. Report I of V: Granular Material Characterization, Department of Civil Engineering, University of Maryland, College Park, 1988.

[11] 董城.路基土動態回彈模量與臨界動應力及路基路面性能研究[D].長沙:中南大學,2014.

[12] TOMMY C, HOPKINS, TONY L, et al. Resilient modulus of Kentucky solis[R]. Lexington: University of Kentucky, 2004.

[13] 劉明. 線性回歸模型的統計檢驗關系辨析[J]. 統計與信息論壇, 2011, 26(4): 21-24.

LIU Ming. Analysis of the relationship among the statistical tests on linear regression model[J]. Statistics & Information Forum,2011, 26(4): 21-24.

[14] 費康,張建偉.ABAQUS在巖土工程中的應用[M].北京:中國水利水電出版社,2010.

[15] 王仕傳.路床加筋力學行為與設計方法的研究[D].上海:同濟大學,2007.