一種新的動態(tài)回彈模量有限元實現(xiàn)方法

寧夏元，董 城，，李志勇，羅衛(wèi)華，梅作舟，冷伍明

（1.湖南省交通科學(xué)研究院，湖南 長沙 410015；2.中南大學(xué) 土木工程學(xué)院，湖南 長沙 410075；3.湖南大學(xué) 土木工程學(xué)院，湖南 長沙 410082；4.西南交通大學(xué) 力學(xué)與工程學(xué)院，四川 成都 610031）

1 引言

目前我國路基路面結(jié)構(gòu)分析是基于靜力學(xué)的彈性層狀體系理論[1－3]，分析采用的參數(shù)是通過靜載試驗獲取，這在低速、低軸載下是合理的。然而，隨著高速、重載公路的迅猛發(fā)展，開展基于動力學(xué)理論的路基路面結(jié)構(gòu)設(shè)計方法成為了必然。Seed等[4]為表征路面回彈變形與路面性能的良好相關(guān)性，提出回彈模量的概念。自此，回彈模量已成為目前國內(nèi)外柔性路面設(shè)計的路基土力學(xué)特征的主要參數(shù)之一[5]。為表征回彈模量與加載應(yīng)力狀態(tài)的相關(guān)性，大量的回彈模量模型被提出[6－10]。根據(jù)所選應(yīng)力變量的不同可分為僅考慮剪切影響的模型[6－7]、僅考慮圍壓影響的模型[8]和綜合考慮剪切和圍壓影響的模型[9－13]。前兩類模型由于僅關(guān)注了剪應(yīng)力或體應(yīng)力單一因素對回彈模量的影響，這與多數(shù)土和粒料回彈模量同時受偏應(yīng)力和體應(yīng)力影響的實際情況不符。因此，在應(yīng)用中存在一定的局限。Uzan[9]在1985年提出了考慮剪應(yīng)力和體應(yīng)力影響的三參數(shù)模型，該模型可以描述粒料的硬化和軟化行為，但存在模量不定值和量綱不一致問題。美國在2004年NCHRP 1-37A項目中提出公路路基設(shè)計的動態(tài)回彈模量模型[10]（后文簡稱為N37A回彈模型）。該模型對Uzan提出的模型進(jìn)行了修正，消除了模量不定值和量綱不一致問題，因而被工程界廣泛接受。后續(xù)許多回彈模量預(yù)估模型均是在N37A回彈模型的基礎(chǔ)上進(jìn)行若干修正[11－13]。

為了利用上述回彈模量模型對路基路面結(jié)構(gòu)進(jìn)行有限元分析，必須采用強(qiáng)健的算法對回彈模量預(yù)估模型進(jìn)行有限元移植。對于N37A回彈模型的有限元實現(xiàn)，已有文獻(xiàn)均采用了基于等效三軸試驗簡化方法[14]，該方法通過等效的軸向應(yīng)力和應(yīng)變推導(dǎo)出等效切線剛度，然后與彈性張量耦合形成等效切線剛度矩陣。該方法推導(dǎo)過程簡單，目前被國內(nèi)外許多學(xué)者所引用[15]，但考慮到該實現(xiàn)方法是基于特殊的等效動三軸試驗方法進(jìn)行的簡化推導(dǎo)，在復(fù)雜應(yīng)力狀態(tài)下會帶來一定的偏差。

因此，需要在廣義線彈性虎克定律的基礎(chǔ)上，嚴(yán)格推導(dǎo)出N37A回彈模型的一致切線剛度矩陣，進(jìn)而編寫用戶自定義材料子程序（UMAT），將回彈模量移植到有限元軟件ABAQUS中。通過簡單軸壓和圍壓加載以及典型路基路面結(jié)構(gòu)變形分析、討論提出的有限元實現(xiàn)方法和已有文獻(xiàn)方法的優(yōu)劣，為基于動態(tài)回彈模量的路基路面結(jié)構(gòu)設(shè)計提供高精度的數(shù)值模擬方法。

2 回彈模量的有限元實現(xiàn)

利用ABAQUS提供的用戶材料子程序（UMAT），可將N37A回彈模型移植到ABAQUS材料庫中。用戶需要通過UMAT定義單元材料積分點(diǎn)的雅可比矩陣，即材料本構(gòu)關(guān)系的切線剛度矩陣。ABAQUS求解時將對每個用戶定義材料單元、材料積分點(diǎn)和平衡迭代都將調(diào)用一次UMAT[16]。

N37A回彈模型的方程為

式中：Mr為回彈模量；θ為體應(yīng)力；τoct為八面體剪應(yīng)力；k1、k2和 k3為路基土應(yīng)力相關(guān)性的回歸系數(shù)；pa為大氣壓力，這里取為100 kPa。

2.1 等效切線剛度矩陣（簡稱為ESTM）

在現(xiàn)有的N37A回彈模型的有限元實現(xiàn)過程中，均采用等效動三軸試驗的簡化有限元實現(xiàn)方法[14]，該方法推導(dǎo)過程簡單。以下對其有限元實現(xiàn)過程做簡要描述。

式（1）中的θ 和 τoct可表示為為偏應(yīng)力張量，其表達(dá)式為I為二階單位矩陣。

式中：S為應(yīng)力張量；

令與（θ，τoct）相對應(yīng)的三軸試驗的軸向應(yīng)力為σa，圍 壓為σr。則θ=σa+2σr，τoct=可求得

在等效三軸試驗應(yīng)力（σa，σr）作用下，相應(yīng)的等效軸向應(yīng)變?yōu)?/p>

式中：σa0和 σr0分別為無輪載作用時路面結(jié)構(gòu)中的垂直應(yīng)力與水平應(yīng)力；ν為泊松比。

定義軸向等效切線剛度 ET為

將式（4）代入式（5）可得

由式（1）求得

在式（7）兩側(cè)同時除以Mr可得

將式（8）代入式（6），最終得到三軸試驗條件下的等效切線剛度 ET計算式為

將該切線模量耦合到彈性張量中，便得到廣義的等效切線剛度矩陣

注意到，式（9）直接過渡到式（10）并無嚴(yán)格證明；對于復(fù)雜的應(yīng)力狀態(tài)，該等效方法是否適用，目前尚未發(fā)現(xiàn)相關(guān)討論。因此，需要推導(dǎo)出UMAT所需的一致性切線剛度矩陣，才能通過對比討論等效切線剛度矩陣帶來的誤差。

2.2 一致性切線剛度矩陣（簡稱為CSTM）

與線彈性本構(gòu)關(guān)系類似，應(yīng)力相關(guān)彈性材料本構(gòu)關(guān)系可以寫成如下形式：

式中： E為應(yīng)變張量；ε為體應(yīng)變，表達(dá)式為ε=tr（E ），令：

其中k=k1(1+ν)，結(jié)合（11）、（12），應(yīng)力相關(guān)彈性本構(gòu)關(guān)系可以簡寫成為

顯然θ、τoct與ξ 之間存在如下函數(shù)關(guān)系：

從而C=C（θ，τoct）可以表示為應(yīng)變的形式：

進(jìn)而式（13）可以表示為

由式（17）、（18），▽CE可以分解成如下形式：

由式（12）可得

由式（16）可得

將式（22）～（25）代入式（21）可得

結(jié)合（14）、（15）和式（20），一致性切線剛度矩陣可以寫成偏應(yīng)力的形式：

利用張量運(yùn)算，式（27）可進(jìn)一步寫為

通過兩種切線剛度矩陣推導(dǎo)過程的對比可知，在等效切線剛度矩陣推導(dǎo)中，式（4）采用Mr來代替從初始應(yīng)力狀態(tài)到當(dāng)前應(yīng)力狀態(tài)的割線模量，使得在初應(yīng)力不為0時得到的等效軸向切線模量為近似值。同時等效切線剛度表達(dá)式（10）中包含了初始應(yīng)力因素，即等效切線剛度矩陣不僅與當(dāng)前應(yīng)力相關(guān)，還與初始應(yīng)力相關(guān)，需要反復(fù)迭代和人為調(diào)整容差來降低迭代次數(shù)，這將給計算過程帶來一定誤差，同時也會帶來收斂性困難。一致切線剛度矩陣由廣義虎克定律推導(dǎo)得到，只與當(dāng)前應(yīng)力相關(guān)，計算更容易收斂。

3 回彈模量數(shù)值實現(xiàn)驗證

3.1 單個單元模擬對比

為驗證N37A回彈模型有限元實現(xiàn)過程中兩種切線剛度矩陣的優(yōu)劣，首先采用單個單元進(jìn)行驗證，單元加載示意圖如圖1所示。模擬中均采用固定載荷步100步，其中軸向載荷σa=100 kPa，圍壓應(yīng)力rσ 分別為1、10、30、50 kPa。

圖1 單元加載示意圖Fig.1 Illustration of element loading

模擬路基土為最佳含水率、壓實度為90%的粉土，N37A回彈模型的參數(shù)如表1所示[17]。

表1 N37A回彈模型參數(shù)Table 1 Parameters of N37A resilient model

圖2給出了不同圍壓下回彈模量及應(yīng)力-應(yīng)變模擬結(jié)果，在不同載荷下兩種切線剛度矩陣計算得到的回彈模量是相同的，與理論解完全吻合；隨著圍壓的增加，CTSM計算得到的應(yīng)力-應(yīng)變曲線始終與解析解吻合，ETSM計算得到的應(yīng)力-應(yīng)變曲線與解析解之間的偏差卻隨著圍壓的增加逐步加大。

圖2 不同圍壓下的模擬結(jié)果Fig.2 Simulated results under different confining pressures

上述結(jié)果可解釋如下：在加載過程中，施加載荷包括軸壓和圍壓，即應(yīng)力狀態(tài)是已知的；而回彈模量的計算僅與當(dāng)前應(yīng)力狀態(tài)有關(guān)。因此，不管用何種切線剛度矩陣均不會影響到回彈模量的計算。然而，應(yīng)變是根據(jù)初始條件通過切線剛度矩陣計算得到，兩種切線剛度矩陣推導(dǎo)過程有很大的差異：ESTM是通過等效切線模量在彈性張量中進(jìn)行推廣得到，這種推廣并沒有得到嚴(yán)格的證明；CSTM基于廣義虎克定律，通過偏應(yīng)力和偏應(yīng)變推導(dǎo)得到的四階張量來進(jìn)行增量求解。ESTM這種等效簡化方法，雖然沒有影響到回彈模量的計算，但計算的軸向應(yīng)力-應(yīng)變曲線與解析解存在一定的偏差。

表2是對上述4種工況應(yīng)力-應(yīng)變曲線最大誤差的統(tǒng)計。由表可見，采用CSTM計算結(jié)果始終與解析解十分吻合，最大誤差為1.95%；，而ESTM在低圍壓（1 kPa）時誤差小于1%，但隨著圍壓的增加，誤差迅速增加；當(dāng)圍壓達(dá)到50 kPa時，誤差增大為13.2%。

表2 有限元解與解析解誤差Table 2 Errors between finite element solution and analytical solution

基于兩種不同的切線剛度矩陣的有限元實現(xiàn)過程中，均采用由應(yīng)變增量與切線剛度矩陣計算應(yīng)力增量，然后再更新應(yīng)力。為了考察兩種切線剛度矩陣對加載步長的依懶性，選取圍壓為50 kPa的加載工況，分析固定載荷步為50、100和1 000步3種情況和自動時間步長進(jìn)行對比分析，結(jié)果如圖3所示。

圖3 不同增量步的模擬結(jié)果Fig.3 Simulation results at different increment steps

由圖3可見，兩種切線剛度矩陣在50、100、1000步工況下，計算結(jié)果對加載步長并不敏感，固定加載步長和自動時間步長的結(jié)果幾乎重合。這表明，在兩種切線剛度矩陣計算得到的應(yīng)力-應(yīng)變曲線對加載步長沒有明顯的依賴性，在實際計算中，采用自動時間步長即可得到滿意的結(jié)果。

上述結(jié)果的原因在于，除了在初應(yīng)力條件下等效切線剛度矩陣需要通過迭代更新應(yīng)力，在每個迭代步中，本構(gòu)關(guān)系仍然是線彈性的。因此，后續(xù)的局部迭代只需要一次即可完成。為驗證上述解釋，將前20個增量步下兩種切線剛度矩陣的等效迭代次數(shù)取出進(jìn)行對比，結(jié)果如圖4所示。

圖4 前20個增量步的等效迭代次數(shù)Fig.4 Equivalent iteration number for the early 20 increment steps

由圖可知，基于兩種切線剛度矩陣的有限元計算所需的迭代次數(shù)大體相同。在增量步初期，等效切線剛度矩陣需要迭代獲取應(yīng)力和應(yīng)變的平衡，隨著初始條件影響的消失，后續(xù)迭代次數(shù)均為1，但一致性切線剛度矩陣能夠獲得更為穩(wěn)定的收斂解。

3.2 實例：路面結(jié)構(gòu)變形分析

上節(jié)討論了單個單元的情形，為進(jìn)一步分析復(fù)雜應(yīng)力狀態(tài)下兩種切線剛度矩陣對計算結(jié)果的影響，以瀝青混凝土路面[18]為例。路面總厚度為69 cm，共分5層，路面頂面作用標(biāo)準(zhǔn)行車荷載[19]。路基路面分析模型示意圖如圖5所示。由于將該問題簡化為平面應(yīng)變問題，標(biāo)準(zhǔn)載荷按照靜力轉(zhuǎn)化的原則進(jìn)行了適當(dāng)轉(zhuǎn)換，轉(zhuǎn)換后輪載為132.3 kPa。計算模型深度取4 m，寬度取8 m，路面結(jié)構(gòu)各層材料的尺寸及參數(shù)如表3所示[18]。計算中引入地應(yīng)力平衡[20]，使自重作用在輪載作用前只影響回彈模量而不影響路基變形。

圖6給出了回彈模量隨深度的變化規(guī)律。兩種切線剛度矩陣計算的結(jié)果十分接近，在路基土的頂端，二者有稍許差異，隨著路基深度的增加，二者趨于重合，這表明，在結(jié)構(gòu)分析中采用兩種切線剛度矩陣所獲得的與應(yīng)力相關(guān)的回彈模量差別較小，亦即采用等效切線剛度矩陣代替一致切線剛度矩陣對回彈模量的影響不大。

圖5 分析模型分層和加載示意圖Fig.5 Illustration of layering and loading of the analyzed model

表3 路面結(jié)構(gòu)材料參數(shù)Table 3 Material parameters of pavement structure

圖6 回彈模量隨路基深度變化曲線Fig.6 Evolution curves of resilient modulus with depth of subgrade

圖7給出了豎向位移隨計算深度的變化曲線。兩種切線剛度矩陣計算得到的路基路面豎向位移在線彈性路面結(jié)構(gòu)層（虛線以上）內(nèi)有3%的偏差，且兩種切線剛度矩陣所獲得的豎向位移隨深度變化曲線幾乎平行。鑒于2種結(jié)構(gòu)分析中除了路基土回彈模量采用不同切線剛度矩陣實現(xiàn)方法外，其他因素都相同。因此，上述偏差是由于采用等效切線剛度矩陣代替一致切線剛度矩陣導(dǎo)致的。由圖7可知，基于一致切線剛度矩陣的豎向位移在相同計算深度均大于基于等效切線剛度矩陣的豎向位移。鑒于材料的剛度越大，提供的變形抗力越大，產(chǎn)生的豎向位移越小，這表明，采用等效切線剛度矩陣代替一致切線剛度矩陣一定程度上增大了路基土變形抗力。

本分析僅考慮了單軸雙輪載荷，對于雙軸和三軸重載情形，由于輪組之間的相互疊加作用，采用等效切線剛度矩陣代替一致切線剛度矩陣可能導(dǎo)致的路基路面響應(yīng)偏差會更大。因此，采用基于當(dāng)前應(yīng)力狀態(tài)的一致切線剛度矩陣有利于提高路基路面結(jié)構(gòu)計算的精度。

圖7 豎向位移隨計算深度變化曲線Fig.7 Evolution curves of vertical displacement with calculated depth

4 結(jié)論

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

（2）單個單元測試結(jié)果表明，等效切線剛度矩陣和一致性切線剛度矩陣對加載子步并不敏感，等效切線剛度矩陣計算誤差隨圍壓的增加而增加，最大誤差在圍壓為50kPa時達(dá)到13.2%；而一致性切線剛度矩陣在相同工況下的最大計算誤差僅為1.95%。

（3）對典型瀝青混凝土路面結(jié)構(gòu)變形的分析結(jié)果表明，等效切線剛度矩陣和一致性切線剛度矩陣對回彈模量的計算完全吻合，但采用一致性切線剛度矩陣計算的路面豎向位移時，在不損失計算效率的前提下，克服了基于等效切線剛度矩陣實現(xiàn)方法在復(fù)雜應(yīng)力狀態(tài)下因模型簡化帶來的偏差，提高了復(fù)雜應(yīng)力狀態(tài)下路面位移計算的精度。

[1]中華人民共和國交通部.JTG D50－2006公路瀝青路面設(shè)計規(guī)范[S].北京:人民交通出版社，2006.

[2]中華人民共和國交通部.JTG D40－2002公路水泥混凝土路面設(shè)計規(guī)范[S].北京:人民交通出版社，2002.

[3]中華人民共和國交通部.JTG D30－2004公路路基設(shè)計規(guī)范[S].北京:人民交通出版社，2004.

[4]SEED H B，CHAN C K，LEE C E.Resilience characteristics of subgrade soils and their relation to fatigue failures in asphalt pavements[C]//Proceedings of International Conference on Structural Design of Asphalt Pavement.Michigan:[s.n.]，1962:611－636.

[5]黃仰賢.路面分析與設(shè)計[M].北京:人民交通出版社，1998:215－220.

[6]THOMPSON M R，ROBNETT Q L.Resilient properties of subgrade soils[J].Journal of Transportation Engineering，ASCE，1979，105(1):71－89.

[7]凌建明，蘇華才，謝華昌，等.路基土動態(tài)回彈模量的試驗研究[J].地下空間與工程學(xué)報，2010，6(5):919－925.LING Jian-ming，SU Hua-cai，XIE Hua-chang，et al.Library research on dynamic resilient modulus of subgrade soil[J].Journal of Underground Space and Engineering，2010，6(5):919－925.

[8]楊樹榮，拱祥生，黃偉慶，等.非飽和粘性路基土回彈模量之研究[J].巖土工程學(xué)報，2006，28(2):225－229.YANG Shu-rong，KUNG Johnson H S，HUANG Wei-hsing，et al.Resilient modulus of unsaturated cohesive subgrade soils[J].Chinese Journal of Geotechnical Engineering，2006，28(2):225－229.

[9]UZAN J.Characterization of granular materials[R]//Transportation Research Record 1022，TRB.Washington D C:National Research Council，1985:52－59.

[10]陳聲凱，凌建明，羅志剛.路基土回彈模量應(yīng)力依賴性分析及預(yù)估模型[J].土木工程學(xué)報，2007，40(6):95－99.CHEN Sheng-kai，LING Jian-ming，LUO Zhi-gang.Stress-dependent characteristics and prediction model of the resilient modulus of subgrade soils[J].China Civil Engineering Journal，2007，40(6):95－99.

[11]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.

[12]NI B，HOPKINS T C，SUN L，et al.Modeling the resilient modulus of soils [C]//Proceedings of the 6th International Conference on the Bearing Capacity of Roads，Railways and Airfields.Lisbon:[s.n.]，2002，(2):1131－1142.

[13]董城，冷伍明，李志勇，等.水泥改良高液限粘土動態(tài)回彈模量試驗研究[J].巖土力學(xué)，2013，34(1):133－138.DONG Cheng，LENG Wu-ming，LI Zhi-yong，et al.Experimental study of dynamic resilient modulus of cement-improved high liquid limit clay[J].Rock and Soil Mechanics，2013，34(1):133－138.

[14]STEVEN W，BARRYarry R，ELI L，et al.Development of design methods for geosynthetic reinforced flexible pavements[R]//FHWA Report Reference:DTFH61-01-X-00068.Bozeman:Cobleigh Hall Montana State University，2004:110－114.

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

[16]徐遠(yuǎn)杰，王觀琪，李健，等.在ABAQUS中開發(fā)實現(xiàn)Duncan-Chang本構(gòu)模型[J].巖土力學(xué)，2004，25(7):1032－1036.XU Yuan-jie，WANG Guan-qi，LI Jian，et al.Development and implementation of Duncan-Chang constitutive model in ABAQUS[J].Rock and Soil Mechanics，2004，25(7):1032－1036.

[17]董城，冷伍明，李志勇.粉土動態(tài)回彈模量試驗研究[J].中南大學(xué)學(xué)報(自然科學(xué)版)，2012，43(12):4834－4839.DONG Cheng，LENG Wu-ming，LI Zhi-yong.Dynamic resilient modulus of silt[J].Journal of Central South University (Science and Technology)，2012，43(12):4834－4839.

[18]胡小弟，孫立軍.瀝青路面結(jié)構(gòu)在非均布荷載作用下的三維有限元分析[J].長安大學(xué)學(xué)報(自然科學(xué)版)，2003，23(3):15－20.HU Xiao-di，SUN Li-jun.Analysis of asphalt pavement structure under non-uniform distributed tire pressure with 3D finite element method[J].Journal of Chang’an University(Natural Science Edition)，2003，23(3):15－20.

[19]陸輝，孫立軍.輪載作用下瀝青路面三維非線性有限元分析[J].土木工程學(xué)報，2004，37(7):64－68.LU Hui，SUN Li-jun.Analysis on asphalt pavement under tyre load by three-dimension finite element method[J].China Civil Engineering Journal，2004，37(7):64－68.

[20]代汝林，李忠芳，王姣.基于 ABAQUS的初始地應(yīng)力平衡方法研究[J].重慶工商大學(xué)學(xué)報(自然科學(xué)版)，2012，29(9):76－81.DAI Ru-lin，LI Zhong-fang，WANG Jiao.Research on initial geo-stress balance method based on ABAQUS[J].Journal of Chongqing Technology and Business University(Natural Science Edition)，2012，29(9):76－81.