考慮非線性接觸及車輛動載鋪裝層動力分析

郝艷廣， 袁龍文， 韓勁龍， 陳金鑫， 陳 璐

(1.中交二航局建筑科技有限公司, 武漢 430000; 2.海工結構新材料及維護加固技術湖北省重點實驗室, 武漢 430040; 3.長沙理工大學土木工程學院, 長沙 410000)

由于公路交通發展迅速，鋼-混組合橋已成為中國公路交通運輸重要組成部分。隨著鋼-混組合橋服役時間延長，鋪裝層與鋼材之間的滑移、開裂日益突出。針對這一現象，有必要進行深入研究，探究鋪裝層與鋼材的滑移、開裂機理。

針對鋪裝層與鋼材之間的滑移、開裂，中外學者已進行大量研究：Stoura等[1]針對大型車-橋系統，提出利用二階常微分方程求解車-橋相互作用問題的動態劃分方法，該方法可減少計算量；Shen等[2]基于模型縮聚理論提出求解車橋耦合問題新理論，并利用MATLAB進行計算分析，驗證了該理論可提高求解精度并減少求解時間；Liu等[3]將移動荷載施加于鋪裝層，計算鋪裝層拉應力，研究重載及溫度對鋪裝層瀝青混合料力學行為的影響；周瑞嬌等[4]利用ANSYS建立車-橋耦合模型，分析車輛速度和路面不平度等因素對橋梁動力響應以及跳車的影響；陳水生等[5]提出采用相位角相干生成車輪相干橋面不平順方法并驗證，結果表明，車輪相干橋面不平順增大橋梁、車輛豎向振動響應，因此有必要考慮各車輪相干性關系。

基于上分析可得，相關研究對車輛動載作用下鋪裝層動力響應研究較少，未能考慮輪胎-鋪裝層非線性接觸，因而存在一些不足之處。利用一定理論建立車輛、橡膠輪胎以及鋪裝層-簡支橋有限元模型，考慮輪胎-鋪裝層非線性接觸以及車輛動載，采用顯式算法求解，通過相關試驗驗證，研究鋪裝層動力響應，并分析路面不平度激勵對鋪裝層動力響應的影響。

1 構建車輛模型

1.1 車輛多體動力學建模

根據車輛振動特點，將其適當簡化。R為慣性參考點，與懸架、車軸固結，如圖1所示。假設車體為剛體B1，利用球較O1約束車體B1與懸架相對自由度，車體可沿OY軸振動；利用旋轉鉸O2約束輪胎與車軸相對自由度，輪胎受動力矩和路面摩擦力作用，將繞OX軸旋轉，使得車輛前進[6-7]。

圖1 車輛示意圖Fig.1 Vehicle diagram

懸架采用線性彈簧阻尼器力元模擬，k表示懸架剛度，c表示懸架阻尼，其作用點為車體質心C1，車體B1受懸架作用。C1初始坐標為q0=(x0,y0,z0)，振動坐標為q=(x,y,z)。懸架對車體B1作用線性彈力F1為

F1=-k(|q|-|q0|)-cq′

(1)

1.1.1 車輛廣義坐標

絕對坐標由剛體質心笛卡爾坐標和繞質心轉動歐拉角構成。選擇笛卡爾坐標系作為R絕對慣性基(O0,e0)，以車體質心C1為基點建立連體基(C1,e1)。車體B1矢徑坐標r1由質心C1相對絕對慣性基(O0,e0)的3個笛卡爾坐標(x1,y1,z1)T表示，車體B1姿態坐標采用相對絕對慣性基(O0,e0)的歐拉角λ1i表示，e0表示基底，其中下標i為歐拉參數編號，i=0,1,2,3。歐拉角坐標陣Λ1為

Λ1=(λ10,λ11,λ12,λ13)T

(2)

車體B1廣義坐標q1為

(3)

1.1.2 無約束車體動力學方程

解除車體B1所有動力學約束，則車體B1無約束動力學方程[7]為

M1q″1=Q1

(4)

(5)

(6)

式中：m1為車體B1質量；E為單位陣；L1為歐拉角坐標陣Λ1引入矩陣；J1為車體B1相對其連體基慣性張量的矩陣；F1和T1為作用于車體B1上全部外力主矢坐標矩陣，車體無約束動力學方程為

Mq″=Q

(7)

1.1.3 車體動力學約束方程

球較O1及旋轉鉸O2構成一組完整約束，其約束方程、速度方程和加速度方程為

Φ1(q,t)=0

(8)

Φ1qq′=-Φ1t

(9)

Φ1qq″=ζ1

(10)

ζ1=-[(Φ1qq′)qq′+2Φ1qtq′+Φ1tt]

(11)

1.1.4 受約束車體動力學方程

對于受約束車體，引入與方程對應的拉格朗日乘子μ。利用高斯加速度變分，可得車體力學微分-代數方程為

(12)

1.2 車輛有限元模型

根據上述分析，基于ABAQUS建立1/4車輛有限元模型。利用連接單元建立懸架與車軸，分別賦予Cartesian屬性和Hinge屬性[8]。車輛有限元模型如圖2所示，車輛參數如表1所示。

表1 車輛參數

圖2 車輛有限元模型Fig.2 Finite element model of vehicle

2 構建輪胎模型

相關研究以質量-彈簧體系模擬橡膠輪胎，與實際存在較大誤差[9-11]。基于ABAQUS建立全鋼絲子午線重載輪胎(輪胎型號11.00R20)，最大氣壓850 kPa，充氣后外緣直徑1 050 mm，充氣后外緣斷面寬282 mm。

2.1 Mooney-Rivlin模型

采用統一數學模型建立橡膠本構關系。假定橡膠不可壓縮并考慮高體積模量，利用應變能密度表征橡膠物理特性[12-14]。由唯相理論可知，Cauchy-Green變形張量不變量I1、I2、I3可表示應變能密度函數，即

(13)

式(13)中：Cij為材料常數，由試驗獲得。

應變能密度函數是主伸長率偶次冪函數，主伸長率與變形張量不變量滿足關系

(14)

(15)

(16)

當橡膠不可壓縮時，I3=1。

應變能密度函數可進一步表示為

W=W(I1,I2,I3,C1,C2,…，Cm,d1,d2,…，dn)

(17)

式(17)中：Ik(k=1,2,3)為應變不變量；Ci(i=1,2,…,m)與dj(j=1,2,…,n)分別為超彈性材料力學特性常數、壓縮特性常數。N=1時，式(13)修改為

(18)

式(18)即為Mooney-Rivlin模型，Jel為橡膠材料變形體積比，此時應變能密度函數為變形張量不變量的一次項函數。Mooney-Rivlin模型降低了單軸拉伸試驗誤差，在小應變范圍擬合效果良好。該模型表達式簡單，參數獲取方便，工程中廣泛采用[15-16]。

2.2 輪胎有限元模型

輪胎主要包括胎面、胎體、帶束層和鋼絲簾線等。為提高計算效率并保持一定精度，對其進行適當簡化：①胎側與其鄰近材料相同的部分合并；②忽略輪胎花紋，減少計算量。輪胎模型及網格劃分如圖3所示。

圖3 輪胎模型及網格劃分Fig.3 Tire model and meshing

橡膠單元為C3D6H和C3D8H，鋼絲簾線單元為SFM3D4R，簾線-橡膠復合材料采用Rebar模型[15-16]。橡膠材料參數及鋼絲簾線材料參數分別如表2、表3所示。表中C10、C01為超彈性材料剪切參數，D1為超彈性材料可壓縮性參數，E為彈性模量，ρ為材料密度，μ為泊松比。

表2 橡膠參數Table 2 Parameters of rubber

表3 鋼絲簾線參數Table 3 Properties of steel cord

3 構建鋪裝層-簡支橋模型

3.1 蠕變模型

瀝青混合料力學性能復雜，與組成材料密切相關。由文獻[15]可知，蠕變試驗可用于測定熱伴瀝青混合料的變形性能，因此采用蠕變模型描述瀝青混合料本構關系。車輛動載作用下，瀝青混合料總應變為

ε(t)=εe+εp+εc

(19)

式(19)中：εe為彈性應變；εp為塑性應變；εc為蠕應變。假定車輛動載作用下瀝青混合料未屈服，則εp=0，式(19)整理為

ε(t)=εe+εc

(20)

溫度恒定(本文中瀝青混合料模擬溫度為20 ℃)，εc為時間t與應力σ的函數，即

(21)

式(21)中：A為冪乘指數，m為等效偏應力階，n為時間硬化模型的總時間階數，均為瀝青混合料蠕變參數，并與溫度相關[15]。

恒載作用下蠕變應變率為

(22)

由式(22)可知，此應變率與時間t有關；非恒載作用下蠕變應變率[15]為

(23)

由式(23)可知，此應變率與時間t無關。車輛動載具有隨機性，故采用與時間t無關的蠕變應變率描述瀝青混合料力學行為。蠕變參數如表4所示[15]。

表4 蠕變參數Table 4 Creep parameters

3.2 鋪裝層-簡支橋有限元模型

以某高速公路鋼-混組合橋為背景，建立鋪裝層-簡支橋有限元模型。鋪裝層尺寸30 m×10 m×0.3 m，如圖4所示；鋪裝層結構如圖5所示；材料參數如表5所示。

表5 材料參數Table 5 Parameters of materials

圖4 鋪裝層-簡支橋模型Fig.4 Model of pavement layer-simply supported bridge

圖5 鋪裝層結構Fig.5 Pavement structure

鋪裝層與鋼材單元均采用C3D8R。鋪裝層沿厚度方向網格劃分如下：上面層為1個網格，網格尺寸為0.025 m；下面層為1個網格，網格尺寸為0.035 m；混凝土層為4個網格，網格尺寸為0.06 m。局部放大的鋪裝層網格劃分圖6所示。

圖6 鋪裝層-簡支橋網格Fig.6 Mesh of pavement layer-simply supported bridge

4 構建車輛-鋪裝層-簡支橋耦合模型

4.1 路面不平度

路面不平度是Gauss平穩隨機過程，一般采用功率譜密度函數表示。利用Fourier逆變換生成路面不平度功率譜密度函數[16]，即

Gd(n)=Gd(n0)(n/n0)-w

(24)

式(24)中：Gd()為位移功率譜密度，m3；n為空間頻率，m-1；n0為參考空間頻率，0.1 m-1；w為擬合指數。GB/T 7031—2005/ISO 8608:1995《機械振動道路路面譜測量數據報告》將路面不平度分為5個級別[17](A～E)。結合工程實際，選取C級路面不平度。

采用余弦疊加法構建路面不平度，即

(25)

式(25)中：r(x)為路面不平度；x為縱向距離；Gd(ni)為功率譜密度；ni為空間頻率；Δn為頻率改變量；θi為[0～2π]內均布的隨機相位角。C級路面不平度如圖7所示。

圖7 C級路面不平度Fig.7 Pavement roughness of class C

4.2 接觸理論

輪胎與鋪裝層接觸狀態復雜，存在較多影響因素。基于Kalkeer三維滾動接觸理論，建立輪胎-鋪裝層接觸關系，其中鋪裝層為主面，胎面為從面，接觸區域網格劃分如圖8所示。

圖8 接觸區域網格劃分Fig.8 Contact area meshing

計算主面與從面接觸單元位移差，以余能最小為控制變量，分別利用法向接觸關系和切向接觸關系求解接觸區域法向力、切向力以及黏著-滑動狀態分布[17-18]，即

(26)

式(26)中：minC1為余能最小值；Ac為可能輪胎-鋪裝層接觸區域；h、μn、pn分別為接觸面法向間隙、法向位移、法向接觸力；wt、μt、μ′t、pt分別為接觸面切向剛性滑移數值、切向位移、前一時刻切向位移、切向接觸力；S為積分面積；g為滑動摩擦系數(與接觸面材料有關)；P為接觸面法向集中力[18]。

假定輪胎-鋪裝層為平面接觸，由以上分析可得，計算時可將任一單元位移與接觸面上作用力表示為

μn=pnAnn,μt=ptAtt

(27)

式(27)中:Ann為在法向單位力產生的法向位移；Att為在切向單位力產生的切向位移。

4.3 耦合振動方程

完成1/4車輛模型與鋪裝層-簡支橋模型的構建，列出二者動力學方程并聯立求解。耦合振動方程[19-20]為

(28)

式(28)中：y、z分別為車輛、橋梁節點位移向量；Fvg為車輛自重等效節點荷載向量；M、C、K為車輛、橋梁節點質量矩陣、阻尼矩陣、剛度矩陣；Fbv、Fvb為輪胎-鋪裝層接觸區域相互作用節點荷載向量，由形函數分配得到。

將車輛模型和鋪裝層-簡支橋模型視為整體，采用直接積分法求解式(28)。直接積分法有顯式和隱式兩種求解方法，本文中采用顯式求解法[21-22]。

5 實例分析

5.1 驗證輪胎模型

利用徑向加載試驗驗證橡膠輪胎有限元模型的可行性。徑向加載試驗如圖9所示，加載-位移曲線如圖10所示。徑向加載試驗釋放輪輞垂向位移，約束其他方向位移；加載速度為0.2 kN/s；氣壓為0.7 MPa。由圖10可知，試驗值與仿真值曲線趨勢接近，最大誤差為4.768%，表明所建橡膠輪胎有限元模型具有一定的可行性[13]。

圖9 徑向加載試驗Fig.9 Radial loading test

圖10 加載-位移曲線Fig.10 Load-displacement curve

5.2 驗證鋪裝層-簡支橋模型

以某高速公路鋼混組合橋為背景，驗證鋪裝層-簡支橋模型。橋梁為簡支橋梁，跨度為30 m。鋪裝層分層鋪筑，由上往下依次為聚合物改性瀝青混凝土GA、環氧瀝青混凝土EA、C50混凝土。利用布設于下面層底部的傳感器測量車輛動載作用下下面層底部縱向應變以及橋梁模態，現場測量如圖11所示，測量用車如圖12所示，應變計規格見表6。

表6 應變計規格Table 6 Strain gauge specification

圖11 現場試驗Fig.11 Field measurement

圖12 試驗用車Fig.12 Test the car

測量步驟：

(1)本次測量共設置4組測量斷面，每組測量斷面間距2 000 mm，其中橫向應變計2組，縱向應變計2組，見圖13。

圖13 傳感器位置Fig.13 Location of sensor

(2)混凝土層：首先將光纖光柵埋入式應變計放置在預先確定的位置，然后澆筑混凝土。

(3)環氧瀝青混凝土層：待混凝土層施工完成后，再將瀝青混凝土埋入式光纖光柵應變計放置在預先確定的位置，然后鋪設環氧瀝青混凝土，最后鋪設聚合物改性瀝青混凝土。

(4)將應變計傳感頭引出線分別連接電源(本次測量采用柴油發電機)和6224動態采集系統相連，信號儀采用INV360智能處理分析儀。開始測量時，試驗車輛預先行駛1～2個回合，待采集器示數穩定后，即可開始正式測量。

(5)本次測量共進行10、15、20、25、30 m/s等6種車輛行駛速度下鋪裝層動力響應的工況。

將下面層跨中底部縱向應變仿真值與測量值進行對比，如圖14所示，誤差分析如表7所示，簡支橋自振頻率仿真值與測量值的對比如表8所示。由圖14、表7和表8可知，仿真值與測量值曲線比較吻合，誤差在容許范圍內，表明所建鋪裝層-簡支橋模型的計算假定和計算結果符合實際情況。

圖14 下面層底部縱向應變Fig.14 Longitudinal strain at the bottom of the lower layer

表7 誤差分析Table 7 The error analysis

表8 自振頻率對比Table 8 Natural frequency contrast

5.3 鋪裝層動力響應分析

5.3.1 無路面不平度

為研究鋪裝層動力響應，選取Vz=30 m/s時，鋪裝層應力響應時程曲線。觀察點為鋪裝層跨中，如圖15所示，各項應力響應如圖16所示。

圖15 觀察點Fig.15 Observation point

由圖16(a)可知，上面層與下面層縱向應力以受壓為主，上面層最大壓應力為0.245 MPa，下面層為0.240 MPa，混凝土層縱向應力以受拉為主，最大拉應力為1.290 MPa。表明混凝土層與鋼材接觸面拉應力明顯。

由圖16(b)可知，瀝青面層與混凝土層橫向應力差別較大。上面層與下面層縱向應力均以受壓為主，上面層最大橫向壓應力為0.277 MPa，下面層為0.267 MPa。混凝土層橫向應力以受拉為主，最大橫向拉應力為0.805 MPa。當車輛接近觀察點時，混凝土層橫向應力出現拉壓交變；車輛遠離觀察點時，無明顯拉壓交變。

由圖16(c)可知，鋪裝層縱向剪應力以拉壓交變為主，各層均表現出相同的變化趨勢。上面層最大拉應力為0.026 MPa，最大壓應力為0.053 MPa；下面層最大拉應力為0.058 MPa，最大壓應力為0.083 MPa；混凝土層最大拉應力0.027 MPa，最大壓應力為0.038 MPa。由上至下，縱向剪應力呈現出先增大后減小的趨勢，下面層出現最大拉應力和最大壓應力。拉應力易導致鋪裝層出現縱向裂縫，在重復交通荷載作用下，鋪裝層縱向剪應力拉壓交變明顯，極易出現疲勞裂縫。

由圖16(d)可知，瀝青層與混凝土層橫向剪應力均表現出拉壓交變，但差別較明顯。上面層最大拉應力為0.008 MPa，最大壓應力為0.006 MPa；下面層最大拉應力為0.015 MPa，最大壓應力為0.010 MPa；混凝土層最大拉應力0.061 MPa，最大壓應力為0.109 MPa。由上至下，橫向剪應力逐漸增大，且混凝土層橫向剪應力明顯大于瀝青面層。當車輛接近觀察點時，混凝土層先壓后拉，瀝青層先拉后壓。

由以上分析可知，混凝土層橫向應力、縱向應力、橫向剪應力受拉明顯，下面層縱向剪應力拉壓交變明顯。

5.3.2 C級路面不平度

為研究路面不平度對鋪裝層動力響應的影響，基于上述分析，選取C級路面不平度下下面層和混凝土層動力響應時程曲線，各項應力時程曲線如圖17所示。

由圖17(a)可知，混凝土層無路面不平度最大縱向拉應力為1.290 MPa，C級路面不平度最大縱向拉應力為2.230 MPa，比無路面不平度增大72.868%，應力時程曲線出現較多尖角，表明振動明顯。

由圖17(b)可知，混凝土層無路面不平度最大橫向拉應力為0.805 MPa，C級路面不平度最大橫向拉應力為2.060 MPa，比無路面不平度增大155.9%漲幅明顯，且C級路面不平度激勵下，橫向應力時程曲線亦存在較多尖角，表現出隨機振動現象。

由圖17(c)可知，混凝土層橫向剪應力出現拉壓交變，無路面不平度最大拉應力0.061 MPa，最大壓應力為0.109 MPa。C級路面不平度最大拉應力為0.067 MPa，最大壓應力為0.231 MPa，分別比無路面不平度增大9.836%，111.927%，壓應力漲幅明顯，且應力時程曲線存在明顯尖角，隨機振動加劇。

6 結論

建立考慮非線性接觸以及車輛動載的橡膠輪胎、鋪裝層-簡支橋有限元模型并驗證模型可行性，得到如下結論。

(1)輪胎徑向加載試驗測量值與仿真值最大誤差為4.768%，加載-位移曲線趨勢接近；下面層底部縱向應變最大誤差為6.112%，應力時程曲線趨勢接近；橋梁自振模態最大誤差為6.724%。表明計算假定和計算方法具有一定的可行性。

(2) 輪胎-鋪裝層非線性接觸以及車輛動載作用下，混凝土層各項應力均表現出明顯拉應力，下面層縱向剪應力拉壓交變明顯，表明鋪裝層易受拉應力影響，從而導致混凝土層與下面層滑移、開裂等病害，且在重復交通荷載作用下，極易出現疲勞損傷。

(3)C級路面不平度激勵下，鋪裝層各項應力時程曲線均表現出非線性增大，除下面層縱向剪應力外，混凝土層各項應力時程曲線均表現出一定的隨機性，表明路面不平度激勵會增大鋪裝層拉應力，增大鋪裝層破壞概率。為保護鋪裝層，施工和運營期應嚴格控制路面不平度。