汪昌勇，馮志強(qiáng)，李九一，陳薈鍵

( 西南交通大學(xué)力學(xué)與工程學(xué)院，成都 610031)

引 言

目前，在數(shù)值模擬方面，對大多數(shù)固體模型進(jìn)行分析的方法通常是有限元法、邊界元法等。在公路工程中，路基檢測的方法一般分為兩種，一種是鉆取樣本檢測方法，另一種是無損檢測法［1］。無損檢測方法中的落錘式彎沉儀( Falling Weight Deflectometer，F(xiàn)WD) 檢測方法，是通過落錘式彎沉儀對路面施加瞬態(tài)載荷( 該載荷可以很好的模擬行車載荷) ，并通過電腦獲取傳感器數(shù)據(jù)，進(jìn)而分析公路質(zhì)量的一種方法［2-3］。落錘式彎沉儀檢測能夠在短時(shí)間內(nèi)獲取數(shù)據(jù)并進(jìn)行分析，適用于對公路的階段性檢測，能夠及時(shí)得出公路質(zhì)量結(jié)果，節(jié)約公路養(yǎng)護(hù)資金，減少交通事故的發(fā)生［4］。高速公路結(jié)構(gòu)一般分為路面層、基層和底基層［5］。按路基性質(zhì)區(qū)分，公路分為剛性路基、柔性路基和半剛性路基［6］。Al －Khoury 等人［7］通過對路基結(jié)構(gòu)的研究，提出了一套新的理論方法即對線彈性多層結(jié)構(gòu)進(jìn)行計(jì)算和分析的譜元法理論。譜元法的基本思路是首先基于連續(xù)介質(zhì)力學(xué)得到結(jié)構(gòu)的受力平衡偏微分方程;然后利用傅里葉變換將時(shí)域內(nèi)的偏微分方程轉(zhuǎn)化為頻域內(nèi)的常微分方程;其次通過引進(jìn)波譜形函數(shù)，把單元的任一位置的位移用結(jié)點(diǎn)位移表示;最后求解常微分方程并代入相應(yīng)的載荷和邊界條件，求出位移以及其它的變量結(jié)果，并將結(jié)果轉(zhuǎn)化為時(shí)域中的結(jié)果［8］。工程中，路基層常常表現(xiàn)出粘彈性，本文以線彈性得到的譜元法理論為基礎(chǔ)，對粘彈性層結(jié)構(gòu)模型進(jìn)行分析，這可為結(jié)構(gòu)的參數(shù)識別方法打下堅(jiān)實(shí)的基礎(chǔ)。

1 譜元法理論

柱坐標(biāo)系下的半空間如圖1 所示，圖中虛線表示理想邊界，該處路表響應(yīng)為零; p( r，t) 為施加的瞬態(tài)載荷，可用點(diǎn)關(guān)于空間和時(shí)間兩個(gè)獨(dú)立函數(shù)的乘積表示，其表達(dá)式為:

對于各向同性線彈性體，以位移表達(dá)的納維方程為:

式中u 為彈性體的位移; ρ 為彈性體的質(zhì)量密度; λ 和μ為拉梅常數(shù)。

圖1 柱坐標(biāo)系下的半空間

根據(jù)斯托克斯－亥姆霍茲矢量分解定理，位移矢量場可用標(biāo)量勢函數(shù)φ( r，t) 和一個(gè)矢量勢函數(shù)珗ψ( r，t) 表示:

1.1 軸對稱結(jié)構(gòu)的求解推導(dǎo)

由于軸對稱性質(zhì)，以勢函數(shù)表示的位移表達(dá)式和應(yīng)力－位移表達(dá)式分別為:

式中u 和w 分別表示水平位移和垂直位移。將式(4) 和式(5) 代入式(2) ，通過傅里葉變換得:

通過求和得到:

對于給定邊界r = R，n = 1，．．．，N; m = 1，．．．，M就足夠能描述模型整體的振動形態(tài)。于是:

Doyle［9］將譜元法運(yùn)用于軸對稱單元，得到了兩種類型的單元，即二節(jié)點(diǎn)層單元和一節(jié)點(diǎn)半空間單元。二節(jié)點(diǎn)軸對稱層單元如圖2 所示。

圖2 二節(jié)點(diǎn)軸對稱層單元

由于二節(jié)點(diǎn)軸對稱層單元具有入射波和反射波的疊加，因此，其位移表達(dá)式為:

僅考慮垂直方向，則有:

式中:

于是可知系數(shù)Amn、Bmn、Cmn和Dmn由上述節(jié)點(diǎn)位移確定，令式(15) 中的4 × 4 的逆矩陣為因于是有:

根據(jù)柯西應(yīng)力原理［10］，正應(yīng)力、剪應(yīng)力與邊界應(yīng)力的關(guān)系為Tk= τkmnm，其中單位矢量n 是垂直于界面且指向外側(cè)的，于是有:

一節(jié)點(diǎn)半空間單元屬于二節(jié)點(diǎn)層單元的特例，如圖3 所示。由于沒有反射波的產(chǎn)生，因此得出:

式中:

式(18) 中的二階矩陣即為一節(jié)點(diǎn)半空間單元的剛度矩陣。

圖3 一節(jié)點(diǎn)半空間軸對稱單元

根據(jù)庫利－圖基的基2 快速傅里葉算法［11］，離散傅里葉變換對為F( t) 和:

式中，k，n = 0，1，．．．，N －1 ; N 是奈奎斯特頻率的采樣數(shù); tk= k·Δt，Δt 是采樣間隔。

載荷分布形態(tài)如圖4 所示。對于一個(gè)圓柱形載荷，它的半徑為a，q = 1，那么載荷的空間分布表達(dá)式可以寫為:

于是根據(jù)傅里葉—貝塞爾理論［12］，可以確定為:

圖4 載荷分布形態(tài)

1.3 剛度矩陣組裝

譜元法剛度矩陣結(jié)構(gòu)單元劃分的示意圖如圖5 所示。

從圖5 可知，譜元法剛度矩陣的結(jié)構(gòu)包含2 個(gè)有限厚度層和1 個(gè)半無限層，分別以2 個(gè)二節(jié)點(diǎn)單元和1 個(gè)一節(jié)點(diǎn)半無限單元組成，一個(gè)單元模擬一整層。譜元法剛度矩陣km，ωn) 的組裝原理類似于有限元法中剛度矩陣的組裝［13］。譜元法剛度矩陣總方程組為:

圖5 結(jié)構(gòu)單元劃分的示意圖

式中:

通過求和并逆變換得:

1.4 譜元法小結(jié)

傳統(tǒng)有限元法對模型劃分單元數(shù)很多，所需計(jì)算時(shí)間較長，而譜元法計(jì)算過程主要在頻域中，避免了遇到無窮積分的計(jì)算難題，且譜元法的二節(jié)點(diǎn)層單元可以模擬整個(gè)路基層，一節(jié)點(diǎn)單元可以模擬半無限路基層，因此在劃分單元上，譜元法優(yōu)于有限元法。

2 粘彈性介質(zhì)

2.1 偏微分算子法

粘性對粘彈性介質(zhì)中波傳播的影響很大［14］。線性粘彈性固體可以利用胡克定律得到應(yīng)力和應(yīng)變關(guān)系，采用偏微分算子法，應(yīng)力和應(yīng)變表達(dá)式為:

將式(26) 傅里葉變換得到:

2.2 伯格斯(Burgers)模型

伯格斯模型是將Maxwell 模型和Kelvin 模型聯(lián)合一起的模型，如圖6 所示。伯格斯模型能很好的模擬粘彈性的特性［15］。

圖6 伯格斯(Burgers)模型

伯格斯模型的應(yīng)力－應(yīng)變表達(dá)式為:

由式(27) 和式(28) 得到:

式中，E*( ω) 為伯格斯模型復(fù)模量。這里假設(shè)材料對體積行為表現(xiàn)為彈性可壓縮( σii= 3Kεii) 以及多維變形的伯格斯行為( Sij= 2μ*( ω) eij) ，于是有:

式中，K 為體積模量。對于粘彈性層剛度矩陣，需要將線彈性公式(17) 中的μ 替換成式(32) 中的μ*( ω) ，然后通過組裝單元剛度矩陣得到其總體剛度矩陣。

2.3 程序計(jì)算流程

程序計(jì)算流程如圖7 所示。

圖7 程序計(jì)算流程圖

首先通過使用編程軟件自編程序代碼;然后通過計(jì)算程序計(jì)算得到模型表面不同位置的垂直方向位移的數(shù)據(jù);最后將數(shù)據(jù)導(dǎo)入繪圖軟件中，獲取數(shù)據(jù)曲線并對其進(jìn)行分析。

3 算 例

Al－Khoury［7］等人通過譜元法對線彈性層的計(jì)算做出了驗(yàn)證，將計(jì)算結(jié)果與有限元軟件CAPA －3D 的計(jì)算結(jié)果作對比，證明了譜元法對線彈性層計(jì)算的精確性和適用性。對于粘彈性層，該方法同樣適用。對粘彈性層的研究，利用伯格斯模型給出了兩個(gè)算例，首先研究泊松比隨時(shí)間變化的情況，即在頻域中泊松比隨頻率變化的情況;其次是研究泊松比為常數(shù)的情況，在這兩種狀況下分別得出瀝青表面位移隨時(shí)間的變化情況。

3.1 算例一

模型分為瀝青層、地基層和底基層，其厚度分別為100 mm、300 mm 和15 000 mm。根據(jù)彎沉儀位移傳感器的位置，分別計(jì)算了模型表面距載荷源中心位置0 mm、300 mm、600 mm、900 mm、1200 mm、1500 mm、1800 mm 的垂直方向位移，其受載荷的時(shí)程曲線如圖8 所示，為50 ms的瞬態(tài)加載，其最大載荷為50 kN，載荷半徑為150 mm。

圖8 載荷時(shí)程曲線圖

圖8 經(jīng)傅里葉正變換得到的頻譜圖如圖9 所示。從圖9 可知，頻率取0 Hz 到150 Hz 就能滿足載荷隨時(shí)間變化的要求。

圖9 載荷頻譜圖

F^m隨m 的變化如圖10 所示。從圖10 可知，隨m的增加，曲線幅值在不斷的衰減，因此取M =1700 已能滿足計(jì)算結(jié)果精確性的要求。為滿足不同頻率波傳播，并使邊界條件R 能充分滿足所有振型情況，在算例中取R=25 m，這就能很好的計(jì)算出瞬態(tài)載荷作用與模型表面的位移。

圖10 傅里葉－貝塞爾系數(shù)分布圖

在載荷時(shí)程曲線中( 圖8) ，取時(shí)間周期T=1 s，載荷樣本點(diǎn)數(shù)為4096，取樣時(shí)間間隔Δt =0.000 244 s，因此取N=4096。模型結(jié)構(gòu)分三層，分別模擬瀝青、基層和底基層，均作為粘彈性材料來研究，其材料參數(shù)見表1。

表1 粘彈性層參數(shù)

這里假設(shè)材料為彈性可壓縮( σii= 3Kεii) 以及變形( Sij= 2u*( ω) eij) 的伯格斯模型。從表1 中瀝青層的數(shù)據(jù)可得到的復(fù)剪切模量以及復(fù)泊松比隨頻率變化的情況，分別如圖11 和圖12 所示。

圖11 復(fù)剪切模量的實(shí)部與虛部圖

從圖11 可知，當(dāng)頻率ω = 0 rad/s 時(shí)，其剪切模量為0，材料處于一種流體的狀態(tài)。從圖12 可知，伯格斯模型呈現(xiàn)出完全不可壓縮的流體狀態(tài)，但隨著頻率的增加，其可壓縮性逐漸增加。

根據(jù)表1 中的參數(shù)數(shù)據(jù)，以及復(fù)剪切模量( 圖11) 和泊松比的變化規(guī)律( 圖12) ，通過程序代碼可計(jì)算出模型表面距載荷源中心不同距離的垂直位移，如圖13 所示。

圖12 復(fù)泊松比的實(shí)部與虛部圖

圖13 模型表面垂直位移時(shí)程曲線圖(泊松比變化)

不同位置最大位移曲線圖如圖14 所示，描述了在不同位置最大位移的趨勢變化，由于軸對稱性質(zhì)，在距載荷源相同的位置其位移相等。

圖14 不同位置最大位移曲線圖(泊松比變化)

3.2 算例二

模型中其他參數(shù)不變，只將瀝青層、地基層和底基層泊松比設(shè)為定值，都為0.45。通過程序計(jì)算得出距載荷源距離分別為0 mm、300 mm、600 mm、900 mm、1200 mm、1500 mm、1800 mm 的瀝青表面垂直位移，其位移時(shí)程曲線和不同位置最大位移曲線分別如圖15、圖16 所示。

從圖15 和16 可知，泊松比設(shè)為定值時(shí)，其位移時(shí)程和不同位置最大位移的變化趨勢皆基本與泊松比變化時(shí)的趨勢一致。

圖15 模型表面垂直位移時(shí)程曲線圖(泊松比為0.45)

圖16 不同位置最大位移曲線圖(泊松比為0.45)

4 結(jié)束語

公路受到高溫時(shí)，路面出現(xiàn)軟化，使得其材料表現(xiàn)出粘彈性。本文以線彈性得到的譜元法理論為基礎(chǔ)，對粘彈性層結(jié)構(gòu)模型進(jìn)行分析。為了便于計(jì)算，將模型結(jié)構(gòu)簡化為三層，分別模擬瀝青層、地基層和底基層，并假定三層均表現(xiàn)為粘彈性，以泊松比隨時(shí)間變化和泊松比不隨時(shí)間變化為例，通過程序計(jì)算得出模型表面距離不同載荷源中心的垂直位移。

由算例一和算例二可知，隨著距載荷源的距離越長，其同一時(shí)刻產(chǎn)生的位移不斷減小，這符合瞬態(tài)力加載時(shí)產(chǎn)生的波在介質(zhì)中傳播不斷衰減的特征; 同時(shí)，在距離載荷源不同位置產(chǎn)生的最大位移時(shí)刻不同，離載荷源越近，產(chǎn)生的最大位移時(shí)刻靠前，離載荷源越遠(yuǎn)，產(chǎn)生的最大位移時(shí)刻靠后，這符合瞬態(tài)動力學(xué)中的滯后現(xiàn)象。由軸對稱性質(zhì)可知，不同位置處的最大位移曲線變化趨勢表明了動力學(xué)問題可以在某一時(shí)刻考慮成為靜力學(xué)問題。通過使用波譜元法對粘彈性模型的位移計(jì)算的研究，可以更好的為路基層的參數(shù)識別提供可靠的理論和計(jì)算上的支撐。