路基壓實(shí)粉質(zhì)黏土動態(tài)回彈模量的四參數(shù)預(yù)估模型及有限元實(shí)現(xiàn)

董 城, 劉文劼, 李 亮

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

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

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

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

1 動態(tài)回彈模量試驗(yàn)及數(shù)據(jù)分析

1.1 試驗(yàn)概況

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

表1 土樣基本物理性質(zhì)指標(biāo)

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

1.2 應(yīng)力狀態(tài)與回彈模量的關(guān)系

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

(a)

(b)

(c)

(d)

2 動態(tài)回彈模量模型改進(jìn)及評價

2.1 已有模型的改進(jìn)及回歸分析

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

(1)

對于不同路基土，剪應(yīng)力對路基土回彈模量的影響程度也是不同的，據(jù)此考慮在N37A模型形式基礎(chǔ)上，引入k4項(xiàng)用于調(diào)節(jié)偏應(yīng)力對回彈模量的貢獻(xiàn)比例，得到的改進(jìn)模型(簡稱N37AP1)具體為

(2)

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

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

表2 預(yù)估模型參數(shù)回歸結(jié)果(k4=0.1)

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

2.2 不同模型預(yù)估結(jié)果對比

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

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

3 N37AP1一致性切線剛度矩陣的推導(dǎo)

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

線彈性本構(gòu)關(guān)系可表示為

(3)

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

類似于線彈性本構(gòu)關(guān)系，非線性彈性本構(gòu)關(guān)系可以利用剛度矩陣Mr寫成

表3 動態(tài)回彈模量預(yù)測與試驗(yàn)結(jié)果對比(k4=0.1)

(4)

令

(5)

式中：k=k1/(1+ν)；k1,k2,k3,k4為N37AP1模型參數(shù)；θ為體應(yīng)力；Pa為標(biāo)準(zhǔn)大氣壓；τoct為八面體剪應(yīng)力。

結(jié)合式(4)和式(5)，動態(tài)回彈模量本構(gòu)關(guān)系可以簡寫為

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

(6)

由式(6)可得體應(yīng)力表達(dá)式

(7)

體應(yīng)力和偏應(yīng)力可以簡寫為

(8)

(9)

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

(10)

顯然θ,τoct與ζ,γ之間存在函數(shù)關(guān)系

(11)

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

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

(12)

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

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

(13)

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

(14)

由定義υ=|ε|，可得

(15)

(16)

(17)

由式(5)可得

(18)

由式(10)可得

(19)

(20)

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

(21)

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

(22)

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

圖2 UMAT實(shí)現(xiàn)流程圖

4 預(yù)估模型有限元實(shí)現(xiàn)驗(yàn)證

4.1 預(yù)估模型的單元驗(yàn)證

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

表4 N37AP1模型的材料參數(shù)

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

圖3 單元加載示意圖

應(yīng)變的解析解可表達(dá)為

(23)

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

圖4 Mr與σ1的關(guān)系曲線

圖5 ε1與σ1關(guān)系曲線

圖6 ε3與σ3關(guān)系曲線

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

4.2 預(yù)估模型的結(jié)構(gòu)驗(yàn)證

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

表5 路基路面結(jié)構(gòu)參數(shù)

圖8為在自重作用下路基路面結(jié)構(gòu)路基土回彈模量的分布云圖，由于路基土采用了應(yīng)力相關(guān)的動態(tài)回彈模量，路基內(nèi)不同位置的回彈模量值各不相同。

圖7 整體網(wǎng)格

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

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

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

(a)豎向分布

(b)橫向分布

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

5 結(jié) 論

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

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

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

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

參 考 文 獻(xiàn)

[1] 羅志剛. 路基與粒料層動態(tài)模量參數(shù)研究[D].上海:同濟(jì)大學(xué), 2007.

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

[3] 曹?，?朱毅,劉云飛,等. 車輛荷載作用下雙層路基層間動應(yīng)力響應(yīng)試驗(yàn)研究[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] 董城,冷伍明,李志勇.粉土動態(tài)回彈模量試驗(yàn)研究[J].中南大學(xué)學(xué)報(bào)(自然科學(xué)版),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] 陳聲凱,凌建明,羅志剛. 路基土回彈模量應(yīng)力依賴性分析及預(yù)估模型[J].土木工程學(xué)報(bào), 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] 凌建明,蘇華才,謝華昌,等. 路基土動態(tài)回彈模量的試驗(yàn)研究[J]. 地下空間與工程學(xué)報(bào),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] 中華人民共和國交通運(yùn)輸部公路路基設(shè)計(jì)規(guī)范: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] 董城.路基土動態(tài)回彈模量與臨界動應(yīng)力及路基路面性能研究[D].長沙:中南大學(xué),2014.

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

[13] 劉明. 線性回歸模型的統(tǒng)計(jì)檢驗(yàn)關(guān)系辨析[J]. 統(tǒng)計(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] 費(fèi)康,張建偉.ABAQUS在巖土工程中的應(yīng)用[M].北京:中國水利水電出版社,2010.

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