基于精細(xì)積分算法的混凝土路面層狀結(jié)構(gòu)動(dòng)力響應(yīng)求解與分析

張宇， 楊林青

(1.河南省收費(fèi)還貸高速公路管理中心航空港管理處, 河南 鄭州 450019；2.廣東技術(shù)師范學(xué)院 天河學(xué)院)

1 前言

隨著中國經(jīng)濟(jì)的不斷發(fā)展，混凝土路面在中國得到越來越廣泛的應(yīng)用。近年來，在中國高等級(jí)公路建設(shè)中，混凝土路面已經(jīng)成為其路面主要結(jié)構(gòu)形式之一。但隨著公路交通量的飛速增長，以及車輛大型化和嚴(yán)重超載等現(xiàn)象的增多，車輛增多、超載和超速等問題開始變得日益嚴(yán)重，越來越引起了人們對(duì)道路安全問題的關(guān)注。路基是路面結(jié)構(gòu)的支撐體，經(jīng)常出現(xiàn)的路面損壞現(xiàn)象大都是路基穩(wěn)定性差、強(qiáng)度不足、在外荷載作用下導(dǎo)致過量變形所致。合理地選擇路基施工措施是獲得堅(jiān)實(shí)而又穩(wěn)定的路基、保證路基路面整體具有良好使用性能的關(guān)鍵。而構(gòu)建瀝青混凝土路面結(jié)構(gòu)在動(dòng)荷載作用下的精確高效數(shù)值模型一直是路面結(jié)構(gòu)體系設(shè)計(jì)參數(shù)反演理論與技術(shù)研究的核心內(nèi)容之一。正演算法還是依據(jù)路面實(shí)測數(shù)據(jù)進(jìn)行反演得到地下結(jié)構(gòu)層設(shè)計(jì)參數(shù)的前提，精確高效的正演模型對(duì)于提高反演精度和效率具有重要意義。

因?yàn)槁访嫠艿暮奢d主要是來自于汽車的運(yùn)動(dòng)荷載，對(duì)道路工程的動(dòng)力響應(yīng)問題已經(jīng)引起人們足夠重視，隨著FWD的廣泛應(yīng)用，為了通過反演計(jì)算得到路面結(jié)構(gòu)的參數(shù)，出現(xiàn)了許多計(jì)算路面彎沉的正演方法，但是大部分計(jì)算方法是把外力當(dāng)成靜荷載來計(jì)算彎沉盆的最大位移。由于用落錘式彎沉儀測定路面的力學(xué)響應(yīng)是在動(dòng)態(tài)荷載脈沖的作用下，并且路面結(jié)構(gòu)，特別是瀝青路面面層在高溫時(shí)呈黏彈性，如果對(duì)這種路面結(jié)構(gòu)仍采用靜荷載作用下彈性層狀理論，計(jì)算結(jié)果必然與實(shí)際結(jié)果有較大誤差。

隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，越來越多的學(xué)者關(guān)注在動(dòng)荷載作用下混凝土路面的動(dòng)力響應(yīng)問題。雖然有限元方法可以模擬復(fù)雜的路面結(jié)構(gòu)，但是其計(jì)算效率和精度嚴(yán)重受網(wǎng)格的影響。Al等開發(fā)了一種精確的、計(jì)算獨(dú)特的分析層狀介質(zhì)動(dòng)態(tài)特性的正演算法，研究了路面在落錘式彎沉儀(FWD)荷載脈沖作用下的響應(yīng)情況；Simon等提出了一種基于譜單元法的動(dòng)態(tài)分析方法來解釋柔性路面的落錘式彎沉儀(FWD)試驗(yàn)。

近年來精細(xì)積分算法被廣泛應(yīng)用到層狀介質(zhì)動(dòng)力響應(yīng)求解中。方宏遠(yuǎn)等采用基于兩端邊值問題的精細(xì)積分方法，分析了電磁波在層狀有耗介質(zhì)中的反射與透射問題，把控制方程按邊值問題求解，可以避免在運(yùn)算過程中出現(xiàn)指數(shù)矩陣相乘的大數(shù)溢出問題，將數(shù)值計(jì)算結(jié)果與廣義傳播矩陣方法加以對(duì)比；Lin等和Han等運(yùn)用精細(xì)積分算法求解得到了層狀地基上任意形狀基礎(chǔ)動(dòng)力響應(yīng)。該文基于傳統(tǒng)譜單元方法，引入精細(xì)積分算法代替?zhèn)鬟f矩陣算法求解路面層狀結(jié)構(gòu)波數(shù)域內(nèi)的動(dòng)力位移響應(yīng)，從而建立路面層狀結(jié)構(gòu)在FWD荷載作用下的動(dòng)力響應(yīng)求解精確算法。通過與已有算法數(shù)值結(jié)果的對(duì)比驗(yàn)證該文算法的準(zhǔn)確性，在此基礎(chǔ)上針對(duì)路面層狀結(jié)構(gòu)的楊氏模量以及厚度進(jìn)行廣泛的參數(shù)分析，以期為工程實(shí)際提供可靠的理論和數(shù)值基礎(chǔ)，為路面結(jié)構(gòu)設(shè)計(jì)和后期的維護(hù)方案選擇提供指導(dǎo)幫助。

2 問題描述

在各種無損檢測方法中，落錘式彎沉儀因其具有類似于卡車車輪施加的高振幅動(dòng)態(tài)載荷的能力，已成為測量路面撓度最廣泛使用的技術(shù)。FWD通過在0.3 m直徑的圓形板上從特定高度扔下一個(gè)落錘，對(duì)路面產(chǎn)生沖擊載荷。撞擊產(chǎn)生的波在路面系統(tǒng)中傳播，由此產(chǎn)生的撓度由n個(gè)檢波器在距離圓形板中心D的地方測量。利用FWD對(duì)瞬態(tài)載荷和撓度進(jìn)行監(jiān)測，測量彎沉盆中心點(diǎn)的位移響應(yīng)。

在路面層狀結(jié)構(gòu)的表面半徑為a的圓形區(qū)域上施加瞬態(tài)荷載Q(t)，采用柱坐標(biāo)系下的軸對(duì)稱幾何形式可以較好地反映FWD測試的力學(xué)幾何特性。該路面由多層具有有限厚度hi的水平層以及底部均質(zhì)半無限空間組成，每層介質(zhì)均為均質(zhì)彈性各向同性屬性。

柱坐標(biāo)系下軸對(duì)稱問題的動(dòng)力平衡方程如下：

(1)

式中：u(r,z,t)和w(r,z,t)分別為徑向和豎向位移；σr(r,z,t)、σθ(r,z,t)和σz(r,z,t)分別為徑向、切向和豎向正應(yīng)力；τrz(r,z,t)為豎直面的切應(yīng)力；ρ為介質(zhì)的質(zhì)量密度。

幾何方程可以表示為：

(2)

式中：εr(r,z,t)、εθ(r,z,t)和εz(r,z,t)分別為徑向、切向和豎向正應(yīng)變；γrz(r,z,t)為豎直面的切應(yīng)變；wθ(r,z,t)為切向轉(zhuǎn)動(dòng)分量。

應(yīng)力-應(yīng)變關(guān)系可以表示為：

(3)

式中：e(r,z,t)=εr(r,z,t)+εθ(r,z,t)+εz(r,z,t)；λ=2vG/(1-2v)，G=0.5E/(1+v)；E、ν和G分別為彈性模量、泊松比和剪切模量。

3 勢能方程

在亥姆霍茲分解中，材料的位移場表示為標(biāo)量勢φ的梯度和矢量勢ψ的旋度之和，如式(4)所示：

u=φ+×ψ

(4)

在軸對(duì)稱運(yùn)動(dòng)中，矢量勢ψ只有一個(gè)分量ψθ。這個(gè)性質(zhì)把問題的解簡化為只對(duì)標(biāo)量勢的解。為了方便標(biāo)記，用ψ代替ψθ。同樣，由于軸對(duì)稱，θ方向的位移分量等于零。用u和w分別表示r和z方向的位移分量，位移分量與勢的關(guān)系由式(5)、(6)表示：

(5)

(6)

將式(5)、(6)代入式(2)，并結(jié)合式(1)、(3)可得：

(7)

(8)

式中：常數(shù)cp=[(λ+2G)/ρ]1/2為壓縮波速；cs=(G/ρ)1/2為剪切波速；式(7)、(8)為兩種基本類型的波：壓縮波(P)和剪切波(S)。

(9)

(10)

(11)

(12)

(13)

設(shè)s=kr，利用鏈?zhǔn)椒▌t，方程(12)可歸納為貝塞爾方程：

(14)

(15)

其中A1為邊界條件確定的常數(shù)，k為徑向波數(shù)。為了使問題離散化，引入徑向邊界條件是必要的。在徑向邊界r=R處(遠(yuǎn)離波源處)，認(rèn)為振蕩振幅消失。因此，式(15)的非平凡解變?yōu)椋?/p>

(16)

這個(gè)條件可以滿足于J0函數(shù)的無窮多個(gè)正根αm。式(16)表示kR=αm，因此k=km=αm/R。m函數(shù)表示為：

(17)

式(17)為波動(dòng)在r方向上的解[方程(12)]，在r=R處消失。每個(gè)函數(shù)對(duì)應(yīng)于第m個(gè)振型。這里可以注意到，r=R處的空間邊界條件不可避免地導(dǎo)致一組離散的正常模態(tài)振動(dòng)波數(shù)。

(18)

其中A2m為常數(shù)，由邊界條件確定，i為復(fù)數(shù)，kpzmn為豎直方向的波數(shù)，表示為：

(19)

式中：cp為壓縮波速；ωn為角頻率。

因此，對(duì)于任意給定的角頻率ωn，式(11)的解為：

(20)

其中Amn為常數(shù)，由邊界條件確定，J0為第零階貝塞爾函數(shù)。

同理可以得到式(10)的解為：

(21)

4 波動(dòng)問題的譜雙重求和

(22)

(23)

N和M的極限值可以由加載脈沖的時(shí)間和空間分布的振幅譜確定，對(duì)N個(gè)頻率的求和可以使用快速傅里葉變換(FFT)，對(duì)M個(gè)波數(shù)的求和可以使用傅里葉-貝塞爾級(jí)數(shù)。傅里葉級(jí)數(shù)的雙求和法是該方法相對(duì)于那些依賴于零到無窮積分的數(shù)值計(jì)算方法的一個(gè)重要的計(jì)算優(yōu)勢。如果系統(tǒng)沒有阻尼或?qū)τ谛∽枘嵊蟹浅＜怃J的峰值，那么這種集成將涉及奇點(diǎn)，并且需要相當(dāng)長的計(jì)算時(shí)間和計(jì)算能力。

在得到φ和ψ的解之后，代入式(5)、(6)中可得軸對(duì)稱介質(zhì)表面位移表達(dá)式為：

(24)

由式(24)可知：M個(gè)波數(shù)之和構(gòu)成波的空間相依性(振蕩形狀)，N個(gè)頻率之和構(gòu)成波的時(shí)間相依性。

5 精細(xì)積分算法

精細(xì)積分是求解一階常微分方程的高精度算法，已被應(yīng)用于結(jié)構(gòu)動(dòng)力分析、控制論、波導(dǎo)等諸多領(lǐng)域。該方法通過將原有計(jì)算層厚度劃分為2N1個(gè)子層，然后再將每個(gè)子層細(xì)分為2N2個(gè)微層，對(duì)于微層，由于其厚度非常小，其層間矩陣可利用冪級(jí)數(shù)展開求解?？梢宰C明級(jí)數(shù)展開的截?cái)嗾`差低于計(jì)算機(jī)的舍入誤差，從而得到幾乎是計(jì)算機(jī)字節(jié)長上精確的數(shù)值解，并利用結(jié)構(gòu)力學(xué)中的多層子結(jié)構(gòu)算法對(duì)這些微段進(jìn)行合并消元，在保證數(shù)值結(jié)果精度的前提下，大大提高運(yùn)算效率。

5.1 區(qū)間矩陣初始化

選取微層的層間矩陣作為初始矩陣進(jìn)行求解，利用如下公式計(jì)算系數(shù)矩陣：

(25)

其中：

(26)

層間矩陣F、G和Q由泰勒級(jí)數(shù)展開計(jì)算得到：

(27)

其中τ為微層的厚度，系數(shù)矩陣計(jì)算公式為：

θ1=-B,γ1=-D,φ1=A

(28)

(29)

(30)

(31)

5.2 區(qū)間矩陣合并

首先對(duì)微層采用如下公式進(jìn)行合并：

(32)

由于微層的區(qū)間矩陣相同，因此可以簡單地對(duì)式(32)執(zhí)行N2次合并循環(huán)即可得到子層的區(qū)間矩陣。

其次在得到子層層間矩陣的基礎(chǔ)上，再采用如下公式對(duì)子層進(jìn)行區(qū)間合并：

(33)

最后，在得到子層區(qū)間矩陣后，進(jìn)一步針對(duì)不同材料屬性的相鄰層逐個(gè)進(jìn)行合并，最終得到均質(zhì)半無限地基之上的所有層的區(qū)間矩陣F、G和Q。

5.3 邊界條件

得到所有層的區(qū)間矩陣F、G和Q之后，其上表面、下表面位移(q0和ql)和應(yīng)力向量(p0和pl)之間的關(guān)系為：

(34)

其中：q0和p0分別為z=0處的位移向量和應(yīng)力向量；ql和pl分別為z=l處的位移向量和應(yīng)力向量。

對(duì)于層狀地基底部(第l+1層)為剛性基巖的情況，需要滿足第l層底部位移為0的邊界條件，即：

ql=0

(35)

由式(34)得到路面位移與應(yīng)力的關(guān)系：

p0=(Q+FG-1FH)q0=E(km)q0

(36)

對(duì)于層狀地基下臥半無限空間的情況，需要考慮無限遠(yuǎn)處的輻射邊界條件。均質(zhì)半無限地基頂部位移和應(yīng)力的關(guān)系為：

pl=R∞ql

(37)

式中：R∞為一個(gè)實(shí)對(duì)稱矩陣，表示半無限空間的剛度矩陣，其計(jì)算可參見文獻(xiàn)[12]。

由式(34)和式(37)對(duì)l層地基底部與半無限地基頂部進(jìn)行組裝，消掉ql和pl可以得到層狀地基表面位移和荷載的關(guān)系為：

p0=(Q+R∞F(GR∞+I)-1FH)q0=E(km)q0

(38)

將式(36)、(38)中的E(km)求逆可得頻率-波數(shù)域內(nèi)層狀地基表面在未知荷載作用下的位移為：

q0=S(km)p0且S(km)=E-1(km)

(39)

式中：S(km)中的元素即為方程(24)中的傳遞函數(shù)。

6 數(shù)值算例

6.1 算法驗(yàn)證

(1)R、M和N0的確定

式(24)中的雙重求和積分涉及到上限的選取，為獲得精確的結(jié)果，選取一路面結(jié)構(gòu)，對(duì)參數(shù)R、M和N0的選取進(jìn)行說明。路面結(jié)構(gòu)包括3層，分別為混凝土面層、基層和底基層，其材料參數(shù)見表1，底部地基假定為均質(zhì)半無限情況，其材料參數(shù)同樣參見表1。

表1 混凝土路面結(jié)構(gòu)和地基材料參數(shù)

利用該文提出的算法，選取如圖1所示的荷載時(shí)程曲線作為FWD荷載，載荷幅值為40 kN，持續(xù)時(shí)間為30 ms，荷載作用于半徑為0.15 m的圓形表面，即落錘的半徑為0.15 m，計(jì)算落錘中心點(diǎn)處的彎沉點(diǎn)位移。

圖1 FWD荷載時(shí)程曲線

首先說明式(16)中R的選取，其中Fourier變換樣本數(shù)為N=2 048個(gè)，波數(shù)樣本數(shù)M=1 500個(gè)，選取不同的R(R=1、10、30、50和60 m)分別進(jìn)行計(jì)算，得到的落錘中心點(diǎn)位移時(shí)程曲線如圖2所示。

圖2 參數(shù)R的確定

由圖2可見：當(dāng)R=50和60 m時(shí)得到的位移時(shí)程曲線基本相同，從而證明了雙重求和積分對(duì)于R有良好的收斂性，在計(jì)算中可以選擇有限的R進(jìn)行計(jì)算。

其次針對(duì)M的選取進(jìn)行說明，選取相同的路面結(jié)構(gòu)模型進(jìn)行計(jì)算，其中Fourier變換樣本數(shù)為N=2 048個(gè)，R=60 m，選取不同數(shù)量的波數(shù)樣本數(shù)M(M=100、500、800、1 000和1 200個(gè))分別進(jìn)行計(jì)算，得到的落錘中心點(diǎn)位移時(shí)程曲線見圖3。

圖3 參數(shù)M的確定

由圖3可見：當(dāng)M=1 000個(gè)和1 200個(gè)時(shí)得到的位移時(shí)程曲線基本相同，從而證明了雙重求和積分對(duì)于M同樣具有良好的收斂性，在計(jì)算中可以選擇有限的波數(shù)樣本數(shù)M進(jìn)行計(jì)算。

最后，分析Fourier變換樣本數(shù)中N0的選取。選取相同的路面結(jié)構(gòu)模型進(jìn)行計(jì)算，其中R=60 m，M=1 500個(gè)，F(xiàn)ourier變換樣本數(shù)為N=2 048。值得注意的是，并不是所有的Fourier變換樣本均需要計(jì)算，選擇合適且有限的前N0個(gè)樣本進(jìn)行計(jì)算，不僅可以得到精確的結(jié)果，而且可以極大地提高求解效率。針對(duì)不同的N0(N0=10、50、80、120和150個(gè))，計(jì)算得到的落錘中心點(diǎn)位移時(shí)程曲線如圖4所示。

圖4 參數(shù)N0的確定

由圖4可見：當(dāng)N0=120個(gè)和150個(gè)時(shí)得到的位移時(shí)程曲線基本相同，從而證明了雙重求和積分對(duì)于N0同樣具有良好的收斂性，在計(jì)算中可以選擇有限的N0進(jìn)行計(jì)算。

(2) 彎沉位移時(shí)程曲線驗(yàn)證

在驗(yàn)證雙重求和積分參數(shù)上限的選取之后，選擇合理的R、M和ωn即可得到高效而精確的結(jié)果。為驗(yàn)證該文算法對(duì)于不同距離點(diǎn)彎沉位移響應(yīng)的準(zhǔn)確性，選擇與上例相同的路面結(jié)構(gòu)和地基進(jìn)行計(jì)算，取R=50 m，km=300，ωn=8 rad/s，計(jì)算距離落錘中心點(diǎn)距離為r=0、0.375、1和1.5 m位置彎沉點(diǎn)的位移響應(yīng)，計(jì)算結(jié)果與文獻(xiàn)[9]進(jìn)行對(duì)比，如圖5所示。

圖5 該文算法結(jié)果與文獻(xiàn)結(jié)果對(duì)比圖

由圖5可見：對(duì)于不同距離情況計(jì)算得到的彎沉位移，兩者基本吻合，從而驗(yàn)證了該文算法的準(zhǔn)確性。

6.2 楊氏模量E對(duì)路面結(jié)構(gòu)動(dòng)力響應(yīng)的影響

路面結(jié)構(gòu)楊氏模量是路面結(jié)構(gòu)設(shè)計(jì)的一個(gè)重要指標(biāo)，其值直接影響著道路的承受能力和使用壽命。為研究其對(duì)路面結(jié)構(gòu)動(dòng)力響應(yīng)的影響，選取3層路面結(jié)構(gòu)進(jìn)行求解分析，包括混凝土面層、基層和底基層，地基假定為均質(zhì)半無限空間，各層材料屬性見表1。

(1) 面層楊氏模量

為分析混凝土面層楊氏模量對(duì)路面動(dòng)力響應(yīng)的影響，將面層的楊氏模量E分別減小50%、增大50%和增大100%進(jìn)行計(jì)算，其他參數(shù)如表1所示。計(jì)算得到的荷載作用點(diǎn)的位移響應(yīng)以及距離荷載作用點(diǎn)2.0 m內(nèi)的位移響應(yīng)最大值如圖6所示。

由圖6(a)可見：針對(duì)不同的楊氏模量E，位移時(shí)程響應(yīng)出現(xiàn)最大位移的時(shí)間基本一致，而最大位移值變化較大，隨著E的增大，最大位移出現(xiàn)明顯的減小現(xiàn)象；由圖6(b)可以看出：面層楊氏模量E的大小，雖然對(duì)荷載施加點(diǎn)的最大位移值影響較大，但對(duì)于距離大于0.5 m的位置影響變小，其影響范圍為0～0.5 m。

(2) 基層楊氏模量

為分析基層楊氏模量對(duì)路面動(dòng)力響應(yīng)的影響，將基層的楊氏模量E分別減小50%、增大50%和增大100%進(jìn)行計(jì)算，其他參數(shù)如表1所示保持不變。計(jì)算得到的荷載作用點(diǎn)的位移響應(yīng)以及距離荷載作用點(diǎn)2.0 m內(nèi)的位移響應(yīng)最大值如圖7所示。

圖6 面層楊氏模量對(duì)路面動(dòng)力響應(yīng)的影響

圖7 基層楊氏模量對(duì)路面動(dòng)力響應(yīng)的影響

由圖7可見：針對(duì)不同的楊氏模量E，其對(duì)路面結(jié)構(gòu)動(dòng)力響應(yīng)的影響與面層基本相似，但是其影響的范圍比面層更大，在距離荷載施加點(diǎn)0～1.0 m范圍內(nèi)基本都有影響，超過1.0 m范圍影響較小。

(3)底基層楊氏模量

現(xiàn)把底基層的楊氏模量E分別減小50%、增大50%和增大100%。利用該文算法，求得路面結(jié)構(gòu)動(dòng)力響應(yīng)荷載作用點(diǎn)處的位移時(shí)程響應(yīng)及距離荷載作用點(diǎn)2.0 m內(nèi)的位移響應(yīng)最大值如圖8所示。

圖8 底基層楊氏模量對(duì)路面動(dòng)力響應(yīng)的影響

由圖8可見：底基層楊氏模量的不同，對(duì)于彎沉位移時(shí)程曲線幅值的影響規(guī)律與面層和基層基本相同，但是其對(duì)最大位移值的影響范圍更大，在0～2.0 m范圍內(nèi)，最大位移曲線基本發(fā)生了平移。

7 結(jié)語

基于傳統(tǒng)的譜單元法，引入精細(xì)積分算法求解頻率-波數(shù)域內(nèi)彎沉點(diǎn)動(dòng)力響應(yīng)，進(jìn)而求解路面層狀結(jié)構(gòu)在FWD荷載作用下的彎沉位移響應(yīng)。該算法采用雙重求和方式，選用有限的積分參數(shù)值代替無限積分，在保證精度的同時(shí)較大地提高了求解效率。通過數(shù)值算例驗(yàn)證了該文算法的準(zhǔn)確性，在此基礎(chǔ)上針對(duì)路面結(jié)構(gòu)層的楊氏模量對(duì)路面彎沉點(diǎn)位移響應(yīng)的影響進(jìn)行了參數(shù)分析。結(jié)果表明：路面結(jié)構(gòu)層的楊氏模量對(duì)彎沉點(diǎn)最大位移值有顯著的影響，隨著楊氏模量的增大，位移最大值明顯減小；不同結(jié)構(gòu)層對(duì)位移最大值的影響范圍不同，面層影響范圍較小，基層影響次之，底基層的影響范圍最大。