牛敏強，鄧明超，蔡 良，苑苗苗

(1.廣東省南粵交通仁博高速管理中心新博管理處，廣東 惠州 516899；2.廣東省南粵交通大豐華高速公路管理中心，廣東 梅州 514329；3.華南理工大學廣州學院，廣東 廣州 501800)

1 前 言

目前，我國現(xiàn)行的公路瀝青路面設(shè)計規(guī)范對路面結(jié)構(gòu)進行力學分析時采用彈性層狀體系理論，該理論將路面材料假定為各向同性體，即在不同的應(yīng)力狀態(tài)下材料的拉伸和壓縮性能是一致的，在材料參數(shù)方面表現(xiàn)為拉壓模量以及泊松比均相同。然而，道路材料拉伸和壓縮力學性能明顯不同，若僅采用單一的壓縮模量表征材料的力學性能將導致路面材料的力學性能被高估，使得路面結(jié)構(gòu)計算結(jié)果與實際工作狀態(tài)嚴重不符，極大影響路面在服役中力學響應(yīng)及長期性能[1-4]。實際上，工程上較多的材料存在拉伸和壓縮性能各異的特征，這常常導致采用該材料修筑的工程結(jié)構(gòu)拉壓力學行為明顯不同[5-8]。

而對于筑路材料，幾乎絕大部分的路用材料也具有類似的性質(zhì)。呂松濤等[9-11]對瀝青混合料和水泥穩(wěn)定碎石的拉壓模量進行測定，研究發(fā)現(xiàn)這兩種材料的拉伸和壓縮模量差別較大，其中壓縮模量約為拉伸模量的1.5～2倍。因此，將道路材料看作拉壓模量不同的材料，在路面結(jié)構(gòu)計算過程中，根據(jù)路面結(jié)構(gòu)中各個積分點的拉壓應(yīng)力狀態(tài)相應(yīng)賦予材料不同的拉壓模量及泊松比，建立符合材料實際服役性能的材料參數(shù)模型，對設(shè)計出合理的路面結(jié)構(gòu)具有重要意義。

基于此，以拉壓不同模量理論為基礎(chǔ)，通過ABAQUS有限元軟件的UMAT用戶子程序功能，將拉壓模量不同的材料本構(gòu)進行二次開發(fā)，通過對比以常規(guī)壓縮模量為材料參數(shù)模型的數(shù)值模擬結(jié)果進行對比，研究拉壓模量不同對路面結(jié)構(gòu)力學響應(yīng)分析結(jié)果的影響。

2 理論基礎(chǔ)

2.1 不同拉壓模量理論基本概念

目前，對于材料拉壓模量不同的研究主要分為兩大類：一種是以復合材料為代表的正交各向異性固體和連續(xù)體；另外一種是對各向同性材料，但拉伸和壓縮式的力學行為存在明顯差異。對于第二種材料本構(gòu)模型，Ambartsumyan[12]進行了系統(tǒng)的研究，并發(fā)展形成了拉壓不同模量的彈性理論。該理論可以采用分段的直線函數(shù)對不同拉壓模量材料的應(yīng)力-應(yīng)變關(guān)系曲線進行近似表達(如圖1所示，其中E+為受拉彈性模量，E-為受壓彈性模量)，該簡化模型具有足夠的精度，完全滿足工程應(yīng)用的要求[13]。本研究也采用該模型進行瀝青路面結(jié)構(gòu)層拉壓模量差異對路面結(jié)構(gòu)力學性能進行研究。

圖1 拉壓不同模量材料的雙直線應(yīng)力-應(yīng)變模型

2.2 基本假設(shè)

(1)研究的對象是連續(xù)的、均勻的彈性體；

(2)在任意應(yīng)力狀態(tài)下只發(fā)生彈性小變形；

(3)各點的計算彈性參數(shù)根據(jù)主應(yīng)力符號的不同，分別賦予不同的彈性模量與泊松比；

(4)假設(shè)μ+/E+=μ-/E-。

3 本構(gòu)關(guān)系與有限元模型的建立

3.1 不同拉壓模量本構(gòu)關(guān)系

分段的直線函數(shù)中，材料的本構(gòu)關(guān)系根據(jù)主應(yīng)力的拉壓不同，彈性模量及泊松比分別取值：當主應(yīng)力為正時(受拉)，彈性模量及泊松比分別為E+及μ+；當主應(yīng)力為負時(受壓)，彈性模量及泊松比分別為E-及μ-。因此，根據(jù)經(jīng)典彈性理論建立以下基于主應(yīng)力方向的二維本構(gòu)方程

(1)

ABAQUS在數(shù)值計算過程中，主程序和用戶子程序UMAT相互之間存在數(shù)據(jù)傳遞的協(xié)同工作過程。其主要工作原理如下：在每個增量步開始時，主程序給節(jié)點施加一個外力增量，然后主程序通過UMAT子程序接口進入UMAT，單元當前高斯積分點相關(guān)的變量的初始值隨之傳遞給UAMT子程序相應(yīng)的變量，ABAQUS會根據(jù)上一增量步提供的切線剛度矩陣計算出節(jié)點位移變量，然后形成新的剛度矩陣并進行迭代計算，直至收斂為止，獲得收斂后的應(yīng)變增量和應(yīng)力增量，最后通過子程序接口將這些變量的更新值返回給主程序，進入下一個增量步進行求解計算。ABAQUS調(diào)用UAMT子程序的詳細過程如圖2所示。

圖2 ABAQUS主程序調(diào)用UMAT流程

以上面層瀝青混合料為例，其二維拉壓不同模量本構(gòu)UMAT的子程序采用Fortran語言編寫主要步驟如下。

(1)定義初始變量。

1DDSDDT(NTENS)，DRPLDE(NTENS)，

STRAN(NTENS)，DSTRAN(NTENS)，

2TIME(2)，PREDEF(1)，DPRED(1)，PROPS(NPROPS)，COORDS(3)，DROT(3，3)，

3DFGRD0(3，3)，DFGRD1(3，3)

REAL(KIND=8)：：ESHE(NTENS)，DESHE(NTENS)，RESHE(NTENS)，VSHE(NTENS)

REAL(KIND=8)：：DSSHE(NTENS)，DSTRES(NTENS)，D(NDI，NDI)

REAL(KIND=8)：：Q11，Q12，Q13，Q21，Q22，Q23，Q31，Q32，Q33，EMOD，ENU，EBULK3，EG2，EG，EG3，ELAM，K1，K2

REAL(KIND=8)：：DDSDDE(3，3)，DSTRESS(3)，EMODcompression，ENUcompression，EMODtension，ENUtension

EMODcompression=PROPS(1)

ENUcompression=PROPS(2)

EMODtension=PROPS(3)

ENUtension=PROPS(4)

(2)定義初始彈性剛度，以層狀彈性體系計算值為初始值

IF(ABS(KINC-1.0).LE.1.0E-3)THEN

EMOD=1500000000.0

ENU=0.3

EBULK3=EMOD/(1-2*ENU)

EG2=EMOD/(1+ENU)

EG=EG2/2

EG3=3*EG

ELAM=(EBULK3-EG2)/3

DO K1=1，NDI

DO K2=1，NDI

DDSDDE(K2，K1)=ELAM

END DO

DDSDDE(K1，K1)=EG2+ELAM

END DO

DO K1=NDI+1，NTENS

DDSDDE(K1，K1)=EG

END DO

ELSE

(3)以X方向主應(yīng)力受壓，Y方向主應(yīng)力受拉為例，求解材料本構(gòu)關(guān)系的雅可比矩陣(DDSDDE矩陣)。

IF(STRESS(1).LT.0..AND.STRESS(2).LT.0.)THEN Q11=EMODcompression/(1.0-ENUcompression*ENUcompression) Q12=EMODcompression*ENUcompression/(1.0-ENUcompression*ENUcompression)

Q13=0.0

Q21=Q12

Q22=EMODcompression/(1.0-ENUcompression*ENUcompression)

Q23=0.0

Q31=0.0

Q32=0.0

Q33=EMODcompression/(1.0+ENUcompression)/2.0

C更新DDSDDE的值

DDSDDE(1，1)=Q11

DDSDDE(1，2)=Q12

DDSDDE(1，3)=Q13

DDSDDE(2，1)=Q21

DDSDDE(2，2)=Q22

DDSDDE(2，3)=Q23

DDSDDE(3，1)=Q31

DDSDDE(3，2)=Q32

DDSDDE(3，3)=Q33

(4)計算應(yīng)力增量DSTRESS，并更新應(yīng)力STRESS。

DSTRESS(1)=DDSDDE(1，1)*DSTRAN(1)+DDSDDE(1，2)*DSTRAN(2)+DDSDDE(1，3)*DSTRAN(3)

DSTRESS(2)=DDSDDE(2，1)*DSTRAN(1)+DDSDDE(2，2)*DSTRAN(2)+DDSDDE(2，3)*DSTRAN(3)

DSTRESS(3)=DDSDDE(3，1)*DSTRAN(1)+DDSDDE(3，2)*DSTRAN(2)+DDSDDE(3，3)*DSTRAN(3)

DO t=1，3

STRESS(t)=STRESS(t)+DSTRESS(t)

END DO

RETURN

END

3.2 有限元計算模型

路面模型在水平方向和深度方向取其有限尺寸，模型的尺寸為5 m(橫向)×5 m(縱向)，各結(jié)構(gòu)層厚度及材料參數(shù)取值見表1。本研究采用拉壓不同模量的模型以及基于壓縮模量的模型進行對比分析。瀝青層、基層、路基應(yīng)用四節(jié)點雙線性平面應(yīng)力四邊形單元(CPS4)進行離散處理，劃分有限元網(wǎng)格。網(wǎng)格劃分的密度選擇自由劃分，模型單元數(shù)量為4 544，節(jié)點數(shù)量為4 680。邊界條件為：對垂直于路線方向兩側(cè)x方向(垂直于路線方向)進行約束；模型底面完全約束。層間接觸條件為完全連續(xù)。為了便于模型計算，輪胎與路面接觸面理想化雙輪垂直均布荷載模型，輪胎半徑為10.65 cm，雙輪中心距31.95 cm，荷載大小為0.7 MPa。

表1 材料參數(shù)及模型尺寸

4 典型瀝青路面結(jié)構(gòu)力學分析

根據(jù)《公路瀝青路面設(shè)計規(guī)范》(JTGD50—2017)，選用如圖3所示的計算點位置進行計算分析。其中A點為車輪中心，B點為車輪跡邊緣，C點為兩個輪胎中心，D點為B點與C點的中點(稱為車輪隙r/4處)。

圖3 力學響應(yīng)計算點位置示意圖

4.1 橫向應(yīng)力分析

選取A、B、C以及D點，分析兩種材料參數(shù)模型下沿道路深度方向的橫向應(yīng)力的變化規(guī)律。4個分析點的橫向應(yīng)力分析結(jié)果如圖4所示。由圖可知，橫向拉應(yīng)力出現(xiàn)在ATB上基層和水泥穩(wěn)定碎石底基層。最大橫向拉應(yīng)力出現(xiàn)在水泥穩(wěn)定碎石底基層層底，且壓縮模量模型的兩個輪胎中心處以下的水泥穩(wěn)定碎石底基層層底最大，其中最大層底拉應(yīng)力為2.757 MPa，拉壓不同模量模型最大層底拉應(yīng)力為2.033 MPa，采用拉壓不同模量模型的最大拉應(yīng)力降低了26.3%。由此可見，基于壓縮模量模型的計算結(jié)果進行路面設(shè)計偏保守。

圖4 橫向應(yīng)力分析

4.2 橫向應(yīng)變分析

選取A、B、C以及D點，分析兩種材料參數(shù)模型下沿道路深度方向的橫向應(yīng)變的變化規(guī)律。4個分析點的橫向應(yīng)變分析結(jié)果如圖5所示。由圖可知， 橫向拉應(yīng)變出現(xiàn)在ATB上基層、 級配碎石下基層、 水泥穩(wěn)定碎石底基層以及土基。最大橫向拉應(yīng)變出現(xiàn)在水泥穩(wěn)定碎石底基層層底，且拉壓不同模量模型的兩個輪胎中心處以下的水泥穩(wěn)定碎石底基層層底最大，其中最大層底拉應(yīng)變?yōu)?08 με，壓縮模量模型最大層底拉應(yīng)變?yōu)?84 με，采用拉壓不同模量模型的最大拉應(yīng)變增加了13%。

圖5 橫向應(yīng)變分析

4.3 剪應(yīng)力分析

選取A、B、C以及D點，分析兩種材料參數(shù)模型下沿道路深度方向的剪應(yīng)力的變化規(guī)律。4個分析點的剪應(yīng)力分析結(jié)果如圖6所示。由圖可知，兩個輪胎中心以下的剪應(yīng)力非常小，最大剪應(yīng)力出現(xiàn)在車輪跡邊緣以下的上面層底，為0.167 8 MPa，此處兩種模型的剪應(yīng)力變化非常小，可忽略不計。

圖6 剪應(yīng)力分析

5 結(jié) 論

采用材料拉壓不同模量理論，結(jié)合ABAQUS有限元UMAT子程序，進行了材料拉壓不同模量模型以及基于壓縮模量模型的瀝青路面結(jié)構(gòu)分析，通過分析路面結(jié)構(gòu)層不同計算點位置的橫向應(yīng)力、橫向應(yīng)變以及剪應(yīng)力，得到以下結(jié)論。

(1)最大橫向拉應(yīng)力出現(xiàn)在水泥穩(wěn)定碎石底基層層底，且壓縮模量模型的兩個輪胎中心處以下的水泥穩(wěn)定碎石底基層層底最大，采用拉壓不同模量模型計算得到的最大橫向拉應(yīng)力降低了26.3%。

(2)最大橫向拉應(yīng)變出現(xiàn)在水泥穩(wěn)定碎石底基層層底，且拉壓不同模量模型的兩個輪胎中心處以下的水泥穩(wěn)定碎石底基層層底最大，采用拉壓不同模量模型計算得到的最大拉應(yīng)變增加了13%。

(3)兩種模型計算得到的不同計算點位置最大剪應(yīng)力變化非常小，分析結(jié)果可忽略不計。