汪昌勇，馮志強，李九一，陳薈鍵

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

引 言

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

1 譜元法理論

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

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

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

圖1 柱坐標系下的半空間

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

1.1 軸對稱結構的求解推導

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

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

通過求和得到:

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

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

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

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

僅考慮垂直方向，則有:

式中:

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

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

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

式中:

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

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

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

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

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

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

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

1.3 剛度矩陣組裝

譜元法剛度矩陣結構單元劃分的示意圖如圖5 所示。

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

圖5 結構單元劃分的示意圖

式中:

通過求和并逆變換得:

1.4 譜元法小結

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

2 粘彈性介質

2.1 偏微分算子法

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

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

2.2 伯格斯(Burgers)模型

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

圖6 伯格斯(Burgers)模型

伯格斯模型的應力－應變表達式為:

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

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

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

2.3 程序計算流程

程序計算流程如圖7 所示。

圖7 程序計算流程圖

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

3 算 例

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

3.1 算例一

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

圖8 載荷時程曲線圖

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

圖9 載荷頻譜圖

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

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

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

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

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

圖11 復剪切模量的實部與虛部圖

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

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

圖12 復泊松比的實部與虛部圖

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

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

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

3.2 算例二

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

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

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

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

4 結束語

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

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