基于PIM的層狀路面結(jié)構(gòu)動力響應(yīng)求解與分析

張海廷，楊林青，郭 芳

(1. 河南水利與環(huán)境職業(yè)學(xué)院，河南 鄭州 450008； 2. 中山大學(xué) 土木工程學(xué)院，廣東 廣州 519082; 3. 廣東技術(shù)師范學(xué)院 天河學(xué)院， 廣東 廣州 510540； 4. 河南農(nóng)業(yè)職業(yè)學(xué)院，河南 中牟 451450)

0 引 言

落錘式彎沉儀(falling weight deflectometer, FWD)是目前應(yīng)用最為廣泛的道路彎沉無損檢測設(shè)備，具有快速、精度高、可重復(fù)性好等優(yōu)點。通過對FWD測得的路面彎沉進行反算分析，可以獲得路面各結(jié)構(gòu)層模量等參數(shù)，從而對路面服役性能評估提供依據(jù)。而最初的路面無損檢測策略主要基于FWD試驗結(jié)果[1-3]，并基于層狀介質(zhì)的靜力計算結(jié)果分析FWD的試驗結(jié)果，從而完全忽略了FWD荷載產(chǎn)生的動力響應(yīng)。因此，建立一種精確的動力響應(yīng)計算方法來分析FWD荷載產(chǎn)生的動力響應(yīng)對路面彎沉分析尤為重要。

開展瀝青混凝土路面結(jié)構(gòu)層模量反算的前提是建立交通荷載作用下路面結(jié)構(gòu)力學(xué)響應(yīng)模型。近年來圍繞FWD荷載作用下路面力學(xué)模型已取得一定成果：豐曉等[4]基于彈性層狀體系理論建立了瀝青路面各層應(yīng)變計算簡式，但沒有考慮瀝青混凝土的黏彈特性；任瑞波等[5]基于傳遞矩陣方法，提出了FWD動荷載作用下多層黏彈性結(jié)構(gòu)動力響應(yīng)的解析解，其中方程的逆變換采用了F.DURBIN的Laplace逆變換的數(shù)值方法[6]；鄧學(xué)鈞等[7]采用四邊形截面軸對稱環(huán)狀單元的有限元法，從一種新的角度重新研究路面層狀結(jié)構(gòu)受水平荷載作用時的剪應(yīng)力問題；B. PICOUX等[8]和A. HAMIM等[9]采用有限單元法建立了一種脈沖荷載作用下柔性路面動力響應(yīng)的數(shù)值模型，其基于有限元法的數(shù)值計算不僅在計算邊界采用人工吸收邊界近似模擬，降低精度，而且計算效率較低；LI Maoyun等[10]采用一種軟計算方法求解得到了FWD動荷載作用下路面彎沉響應(yīng)；S. GRENIER等[11]和R. AL-KHOURY等[12]利用譜單元法建立了FWD動荷載作用下柔性路面結(jié)構(gòu)響應(yīng)模型，其頻率波數(shù)域內(nèi)二階常微分方程的解采用了一種近似于傳遞矩陣方法的公式計算；黃兵等[13]和宋金華等[14]建立有限元模型分析了半剛性基層瀝青路面結(jié)構(gòu)動力響應(yīng)；黃志義等[15]采用Witczak 黏彈性模型動態(tài)剪切模量建立了非均布移動荷載作用下瀝青路面結(jié)構(gòu)瞬態(tài)動力分析的三維有限層模型。近年來很多關(guān)于層狀地基的研究表明，地基的非均質(zhì)特性對地基動力響應(yīng)有顯著的影響[16-17]，在實際工程中應(yīng)該予以考慮。

精細積分方法是一種求解一階偏微分方程的高精度數(shù)值算法，尤其適用于計算層狀結(jié)構(gòu)動力響應(yīng)，其數(shù)值精度可達計算機意義的解析解，目前已經(jīng)應(yīng)用于層狀無限地基地震響應(yīng)、無阻尼層狀地基動力響應(yīng)計算[18-19]和層狀介質(zhì)探地雷達電磁波反射信號等領(lǐng)域[20]。王妍等[21]基于精細積分算法和混合變量算法提出了三維路面層狀結(jié)構(gòu)動力響應(yīng)的求解方法。但該方法采用數(shù)值積分變換，雖然用有限積分代替了無窮積分，其求解效率依然較低。筆者應(yīng)用譜單元法，引入精細積分方法代替?zhèn)鹘y(tǒng)譜單元法中的兩節(jié)點單元和單節(jié)點單元，并采用Burgers模型模擬瀝青混凝土材料的黏彈特性，建立了一種交通荷載作用下路面動力響應(yīng)計算模型，計算FWD動荷載作用下瀝青混凝土路面彎沉響應(yīng)，并將結(jié)果與彈性模型以及實測結(jié)果進行對比，驗證模型的精度以及對實際工程的適用性。

1 軸對稱體系波動方程

FWD對路面施加動荷載的數(shù)值模擬可以作為一個軸對稱問題進行求解分析，路面坐落在層狀半無限路基上，如圖1。軸對稱豎直向荷載施加在路面直徑為d的圓盤區(qū)間內(nèi)(圖2)，其表達式為:

X(r,t)=Y(r)Z(t)

(1)

式中：S(r)為軸對稱豎直向荷載空間分布；Z(t)為軸對稱豎直向荷載隨時間的變化。采集點距離圓盤中心的距離用x表示，建立坐標系(圖1)，坐標原點置于圓盤中心。此時各向同性線彈性材料的運動方程可以表示為：

(2)

式中：向量u為位移函數(shù)；ρ為質(zhì)量密度；?為向量微分算子；?·u為u的散度；?2u為u的拉普拉斯算子向量；λ和μ是拉梅常數(shù)，其與楊氏模量E和泊松比ν的關(guān)系為：

(3)

圖1 落錘彎沉儀系統(tǒng)示意Fig. 1 Schematic diagram of the falling weight deflectometer system

圖2 軸對稱豎直向荷載Fig. 2 The axisymmetric vertical load

在Helmholtz分解中，地基的位移場可以表示為標量勢φ的梯度和向量勢ψ的旋度之和，即：

u=?·φ+?×ψ

(4)

在軸對稱運動中，矢量勢ψ只有一個分量ψθ。同樣由于是軸對稱問題，在θ方向上的位移分量等于0。用ur和uz分別表示r和z方向上的位移分量，位移分量與勢的關(guān)系為：

(5)

(6)

引入相關(guān)的應(yīng)力-位移關(guān)系可得軸對稱波動方程：

(7)

(8)

式中：常量cP和cS分別為縱波波速和橫波波速，其表達式如式(9)、式(10)：

(9)

(10)

式(7)、式(8)表示兩種基本類型的波：壓縮波(P波)和剪切波(S波)，兩者相互獨立，可由其波速cP或cS表示其特征。

2 勢能方程

勢函數(shù)φ和ψθ可以從基本運動方程(7)～(8)利用分離變量方法進行求解。通過傅里葉變換將基本運動方程(7)～(8)從時域轉(zhuǎn)化到頻域：

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

式中:ωn為角頻率,n=1,2,…,N。

3 波動的譜雙重求和分析

(19)

(20)

從加載脈沖的時間和空間分布的幅值譜可以確定N和M，對N頻率的求和可以通過快速傅里葉變換完成，而對M波數(shù)的求和可以通過傅里葉-貝塞爾級數(shù)完成。傅里葉級數(shù)的二重求和方法構(gòu)成了該方法相對于那些依賴于零和無窮之間積分的數(shù)值計算方法的基本計算優(yōu)勢。如果系統(tǒng)沒有阻尼，或者小阻尼的峰值非常尖銳，那么這種集成涉及到奇異性，并且它需要相當多的計算時間和容量。

將通解〔式(19)、式(20)〕代入式(6)，然后考慮邊界條件，軸對稱路面的豎直向位移可以表示為：

(21)

1)計算屬性矩陣：利用拉梅常數(shù)λ和μ以及波數(shù)km，通過式(22)計算當前層的屬性矩陣K11、K12、K21和K22:

(22)

2)計算特征矩陣：利用計算得到的屬性矩陣，通過式(23)求解當前層的特征矩陣:

(23)

3)將當前層的厚度劃分為2N1(N1為任意正整數(shù))個相同厚度的子層，則每個子層的厚度為η=hX/2N1。然后，將每一個子層進一步劃分為2N2(N2亦為任意正整數(shù))個微細薄層，則每一個微細薄層的厚度為τ=η/2N2。由此可見，對于微細薄層，其厚度τ非常小，所以與微細薄層對應(yīng)的系數(shù)矩陣ML(τ)、MG(τ)和MQ(τ)可由Taylor級數(shù)展開近似計算。

(24)

經(jīng)過N2次的合并可以得到子層的系數(shù)矩陣。

5)子層和地基層合并：利用計算得到的子層系數(shù)矩陣，通過式(25)進一步合并得到整個層狀地基的系數(shù)矩陣:

(25)

6)底部邊界條件：針對層狀地基底部半無限空間情況，利用求解得到的半無限地基特征矩陣Δ，對其求解特征值問題：

ΔΦ=ΦΛ

(26)

式中：Λ為特征值矩陣；Φ為特征向量矩陣，二者寫成分塊矩陣的形式為：

(27)

(28)

式中:λj的實部均為正值，j=1,2,則半無限地基表面的邊界條件為:

pl=K∞ul

(29)

為了便于計算，對于實際應(yīng)用的豎直向荷載X(r,t)=Y(r)Z(t)，其產(chǎn)生的豎直向位移可以通過式(21)轉(zhuǎn)化為式(30)進行求解：

(30)

(31)

式中：k和n的取值范圍為(0,…,N-1)；N為FFT的頻率樣本數(shù)；Δt為采樣間隔，且tk=k·Δt。

(32)

式中：R為結(jié)構(gòu)的徑向距離(R相對較大，以確保uz在方程中等于零)。

通過傅里葉-貝寒爾級數(shù)可表達出Y(r)，如式(33)：

(33)

4 黏彈性模型

經(jīng)室內(nèi)試驗驗證，各種黏彈性模型中，Burgers模型能較真實地反映瀝青混合料的力學(xué)性能，故筆者采用如圖3的Burgers黏彈模型。此模型集合了Kelvin和Maxwell兩種模型，Maxwell模型由一個線彈性常數(shù)為E1的胡克單元和一個黏性系數(shù)為η1的牛頓單元串聯(lián)組成；Kelvin模型由兩個常數(shù)分別為E2和η2的胡克單元和牛頓單元并聯(lián)組成。

圖3 Burgers黏彈模型Fig. 3 Burgers viscoelastic model

Burgers模型在單軸壓縮或拉伸試驗下的應(yīng)力-應(yīng)變關(guān)系為：

(34)

對式(34)進行傅里葉變換，則有：

(35)

通過類比關(guān)系可得Burgers黏彈模型復(fù)模量為：

E*(ω)=

(36)

該復(fù)模量的值可以通過圓柱體樣本的單軸掃頻實驗獲得，也可通過結(jié)構(gòu)動力無損檢測獲得。將此復(fù)模量替換式(3)中的E即可獲得復(fù)拉梅常數(shù)。

利用式(36)可以得到復(fù)拉梅常數(shù)λ*(ω)和μ*(ω)，如式(37)、式(38)：

λ*(ω)=K-2μ*(ω)/3

(37)

(38)

式中：K為體積模量。將式(37)和式(38)中的λ*(ω)和μ*(ω)替換式(22)中的λ和μ即可。

5 數(shù)值算例

5.1 頻域內(nèi)計算結(jié)果驗證

為驗證該文中算法頻域內(nèi)結(jié)果的準確性，選擇半無限空間的單層地基進行分析。單層地基和半無限地基的剪切波速和質(zhì)量密度之比分別為0.5和1.0，泊松比均為0.33。在地基表面半徑為0.070 5h的圓盤區(qū)間施加豎直向簡諧荷載，計算得到ωnh/cS=7.0(cS表示單層地基橫波波速)時豎直向位移沿r的分布，見圖4。

圖4 豎直向位移沿r的分布Fig. 4 Distribution of the vertical displacement along r direction

圖4中，將計算結(jié)果與文獻[23]進行對比，兩者完全擬合，從而驗證了文中算法頻域計算結(jié)果的準確性。值得注意的是，文獻[23]中算法在求解此類問題時需要進行無限域的積分，為了足夠的精度，需要完成相當大的計算量，而文中算法有較高的求解效率，無需求解無窮積分，在獲取較高精度的同時僅需較小的計算量。

5.2 時域內(nèi)計算結(jié)果驗證

為進一步驗證文中算法在時域內(nèi)計算結(jié)果的準確性，選取一高模量瀝青混凝土路面結(jié)構(gòu)進行計算。該路面結(jié)構(gòu)層共包括3層(如圖5)，即高模量瀝青混凝土面層(HMAC層)、基層和底基層，其材料屬性(彈性模量E、泊松比ν、密度ρ和阻尼比ξ)和厚度見表1。為與文獻[11]的結(jié)果進行對比，對于路面結(jié)構(gòu)的黏彈屬性不采用Burgers模型進行模擬，而是通過式(39)對復(fù)模量進行計算：

E*=E(1+2iξ)

(39)

表1 瀝青混凝土路面結(jié)構(gòu)材料屬性Table 1 Properties of the asphalt concrete pavement structure

圖5 瀝青混凝土路面結(jié)構(gòu)層(單位：cm)Fig. 5 The asphalt concrete pavement structure layer

施加的FWD荷載幅值為40 kN，時長為30 ms，荷載時程曲線見圖6，荷載施加在直徑為30 cm的圓盤區(qū)域，位移接收點距離該圓盤區(qū)域圓心的距離分別為0、25、37.5、50、62.5、75、100、125、150、200 cm。傅里葉變換樣本數(shù)N=2 048，波數(shù)樣本數(shù)M=1 500，式(33)中的R=150 m。通過文中算法計算得到各個彎沉點的位移時程曲線如圖7，并與文獻[11]結(jié)果進行對比。由圖7可見，兩者基本完全擬合，從而驗證了文中算法的準確性。此外，文獻[11]只給出了距離為0、37.5、100、150 cm彎沉點的時程曲線，筆者將其它各彎沉點的時程曲線亦全部繪出。

圖6 FWD荷載時程曲線Fig. 6 The time-history curve of FWD load

圖7 FWD動荷載作用下各彎沉點位移時程曲線Fig. 7 The time-history curve of every deflection point displacement under FWD dynamic load

5.3 黏彈解和彈性解路表彎沉比較

選擇三層瀝青路面結(jié)構(gòu)，如圖8，瀝青混凝土路面層厚10.0 cm，材料黏彈性性質(zhì)采用Burgers模型來描述，其中E1=1 260 MPa，E2=630 MPa，η1=31.1 MPa，η2=23.3 MPa，密度ρ=2 400 kg/m3，泊松比ν=0.35?；鶎雍吐坊鶠閺椥泽w，基層厚度為30.0 cm，路基為半無限空間，兩者的彈性模量分別為1 090、196 MPa，泊松比分別為0.25和0.40，密度分別為2 300、1 900 kg/m3。采用落錘彎沉儀對同一路基進行實驗，采集點距離落錘點的距離分別為0、30、60、90、120、150、180 cm，落錘的半徑為15 cm。

圖8 瀝青路面結(jié)構(gòu)層(單位:cm)Fig. 8 The asphalt pavement structure layer

在FWD動荷載作用下(圖9)，利用文中算法對路表彎沉進行求解。時長為0.06 s，傅里葉變換樣本數(shù)為N=2 048，波數(shù)樣本數(shù)M=1 500，式(33)中的R=150 m。

圖9 FWD試驗荷載時程曲線Fig. 9 The time-history curve of FWD test load

為了分析文中算法黏彈模型的優(yōu)越性，將瀝青混凝土材料屬性退化為線彈性后同樣利用文中算法進行求解，瀝青混凝土路面層彈性模量選為E1，將兩者的計算結(jié)果(r=0處)與落錘彎沉儀測得的實際位移時程曲線進行對比，如圖10。由圖10可以看出，黏彈解相比彈性解與實驗結(jié)果擬合得更好，尤其是曲線的下降段，由于黏彈屬性，路面位移回落的時間相比彈性解滯后，這也體現(xiàn)了黏彈材料本身的一個特點，也更加符合路面材料的實際情況。所有采集點處的最大位移曲線如圖11。由圖11可以看出，黏彈解得到的最大位移均略大于彈性解，更接近于實驗值。研究結(jié)果驗證了筆者提出的黏彈模型在求解瀝青混凝土路面問題的準確性，用此模型來模擬此類路面更符合實際情況。

圖10 FWD動荷載作用下落錘中心(r=0)彎沉位移時程曲線Fig. 10 The time-history curve of deflection basin displacement at r=0 under the FWD dynamic load

圖11 FWD動荷載作用下采集點彎沉最大值Fig. 11 The maximum deflection value of the acquisition point under the FWD dynamic load

5.4 層厚對路面位移響應(yīng)的影響

應(yīng)用筆者提出的計算模型，針對瀝青層和路基層厚度對路面位移響應(yīng)的影響進行參數(shù)分析。選擇三層路面結(jié)構(gòu)層為例，瀝青混凝土路面層為黏彈屬性，其黏彈參數(shù)分別為E1=8.1 GPa、E2=4.05 GPa、η1=0.2 GPa、η2=0.15 GPa、K=8.1 GPa，密度ρ=2 300 kg/m3。路面基層和路基為彈性層，彈性模量分別為10.8 、0.1 GPa，泊松比分別為0.35和0.40，密度分別為2 100 、1 800 kg/m3。

為分析瀝青混凝土路面層厚度對路面位移響應(yīng)的影響，設(shè)定路面基層的厚度為0.6 m，路基為半無限地基，瀝青混凝土路面層厚度分別為5、10、15 cm，在路面施加的FWD荷載幅值為40 kN，時長為30 ms，荷載時程曲線如圖6。荷載施加在直徑為30 cm的圓盤區(qū)域，位移接收點距離該圓盤區(qū)域圓心的距離分別為0、25、37.5、50、62.5、75、100、125、150、200 cm。傅里葉變換樣本數(shù)N=2 048，波數(shù)樣本數(shù)M=1 500，式(33)中的R=150 m。計算得到的彎沉盆中心位移幅值曲線如圖12。

圖12 瀝青混凝土路面層厚度對彎沉盆位移最大值影響Fig. 12 The influence of thickness of asphalt concrete layer on the maximum value of the deflection basin displacement

由圖12可見，瀝青混凝土路面層厚度對彎沉盆中心的位移幅值有顯著的影響。隨著厚度的增大，彎沉點處的位移均減小，但在荷載施加位置的位移減小幅度最?。浑S著距離增大，瀝青混凝土路面層厚度對位移幅值的影響呈先變大后逐漸減小趨勢。

計算得到的彎沉盆時程曲線如圖13。由圖13可見，瀝青混凝土路面層厚度對位移時程的影響主要集中在幅值點附近，在位移的上升段和衰減段的影響均比較小，而其對彎沉點達到位移最大值的時刻基本沒有影響。

圖13 瀝青混凝土路面層厚度對彎沉盆位移時程影響Fig. 13 The influence of thickness of asphalt concrete layer on the time history of the deflection basin displacement

為分析路面基層對路面位移響應(yīng)的影響，設(shè)定瀝青混凝土路面層厚度為0.1 m，路基為半無限地基，路面基層的厚度分別為0.3、0.6、0.9 m，施加相同的FWD荷載，計算得到的彎沉盆中心位移幅值曲線如圖14。由圖14可見，路面基層厚度對彎沉點位移幅值的影響比瀝青混凝土路面層厚度對其的影響更加明顯，其在荷載施加位置的影響最大。隨著路面基層厚度的增加，位移幅值顯著減??；隨著距離的增加，影響逐漸減小。

圖14 路面基層厚度對彎沉盆位移最大值影響Fig. 14 The influence of thickness of subbase layer on the maximum value of the deflection basin displacement

計算得到的彎沉盆位移時程曲線如圖15。由圖15可見，路面基層厚度對彎沉點位移時程曲線的影響規(guī)律與瀝青混凝土路面層厚度的影響規(guī)律基本相同，只是影響幅度更大一些。

圖15 路面基層厚度對彎沉盆位移時程影響Fig. 15 The influence of thickness of subbase layer on the time history of the deflection basin displacement

6 結(jié) 論

筆者基于精細積分算法，提出了一種求解FWD荷載作用下路面動力響應(yīng)的數(shù)值模型，對于瀝青混凝土路面，采用Burgers模型進行模擬，考慮了瀝青混凝土路面的黏彈特性。

頻域和時域數(shù)值算例均驗證了文中算法的正確性，同時將文中算法的黏彈性解與彈性解和實測結(jié)果進行對比。結(jié)果表明：與彈性模型解相比，黏彈模型解計算得到的動力響應(yīng)時程曲線與實測結(jié)果更加擬合。黏彈模型計算得出的最大位移峰值曲線與實測結(jié)果基本擬合，而彈性解計算得出的結(jié)果均小于黏彈解和實測值。研究結(jié)果一方面驗證了筆者建立路面層狀結(jié)構(gòu)動力響應(yīng)計算模型的精度以及對實際工程的適用性，另一方面也驗證了在進行路面結(jié)構(gòu)動力分析時應(yīng)予以考慮瀝青混凝土路面結(jié)構(gòu)的黏彈特性。

利用筆者提出的算法，針對瀝青混凝土路面層和路面基層厚度，對路面動荷載作用下的位移進行廣泛的參數(shù)分析。結(jié)果表明，路面結(jié)構(gòu)層的厚度對于動力響應(yīng)的影響均比較明顯，尤其是路面基層厚度影響最明顯。