材料拉壓彈性模量差異對路面倒裝結構力學性能的影響

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

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

1 前 言

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

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

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

2 理論基礎

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

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

圖1 拉壓不同模量材料的雙直線應力-應變模型

2.2 基本假設

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

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

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

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

3 本構關系與有限元模型的建立

3.1 不同拉壓模量本構關系

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

(1)

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

圖2 ABAQUS主程序調用UMAT流程

以上面層瀝青混合料為例，其二維拉壓不同模量本構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方向主應力受拉為例，求解材料本構關系的雅可比矩陣(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)計算應力增量DSTRESS，并更新應力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(縱向)，各結構層厚度及材料參數(shù)取值見表1。本研究采用拉壓不同模量的模型以及基于壓縮模量的模型進行對比分析。瀝青層、基層、路基應用四節(jié)點雙線性平面應力四邊形單元(CPS4)進行離散處理，劃分有限元網(wǎng)格。網(wǎng)格劃分的密度選擇自由劃分，模型單元數(shù)量為4 544，節(jié)點數(shù)量為4 680。邊界條件為：對垂直于路線方向兩側x方向(垂直于路線方向)進行約束；模型底面完全約束。層間接觸條件為完全連續(xù)。為了便于模型計算，輪胎與路面接觸面理想化雙輪垂直均布荷載模型，輪胎半徑為10.65 cm，雙輪中心距31.95 cm，荷載大小為0.7 MPa。

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

4 典型瀝青路面結構力學分析

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

圖3 力學響應計算點位置示意圖

4.1 橫向應力分析

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

圖4 橫向應力分析

4.2 橫向應變分析

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

圖5 橫向應變分析

4.3 剪應力分析

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

圖6 剪應力分析

5 結 論

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

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

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

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