分層巖土介質(zhì)振動響應的傳遞矩陣法求解

畢戰(zhàn)歌， 劉 軍，2

（1.河海大學土木與交通學院安全與防災工程研究所，南京 210098；2.河海大學巖土力學與堤壩工程教育部重點實驗室，南京 210098）

分層巖土介質(zhì)的動態(tài)力學特性研究一直是國內(nèi)外學者關(guān)注的熱點之一. 在巖土工程、防護工程和爆破工程中，分層巖土介質(zhì)彈性半空間的地表振動響應研究具有重要工程意義. 許多學者在分層介質(zhì)的動力奇異解問題上做了研究，F(xiàn)ranssens［1］導出了多層半無限介質(zhì)的平面問題的時域奇異解；Apsells［2］和Kauselle［3］分別導出了多層半空間問題的頻域奇異解；孫成禹和張立［4］求解了與水平面呈任意夾角的表面線源作用于彈性半空間的拉普拉斯-傅里葉雙積分變換解；奚亞男［5］在勢函數(shù)理論基礎(chǔ)上，運用Laplace變換、傅里葉-貝塞爾積分等數(shù)學手段，得出了地基表面各點位移的形式解. 許多學者對層狀介質(zhì)的振動響應問題從不同的角度上進行了研究. Piotr和Neves［6］等提出了一種基于小波變換的方法來分析和優(yōu)化多層介質(zhì)在移動荷載作用下動力響應的方法. 在試驗研究方面，Wang［7］等利用SHPB試驗裝置，建立了一種適用于層狀介質(zhì)動態(tài)性能研究和數(shù)據(jù)分析的試驗方法. 在數(shù)值模擬方面，Kumar［8］等利用考慮投射邊界的有限元法，對土-巖和巖-巖地基系統(tǒng)在豎向振動模式下對塊體基礎(chǔ)動力響應的影響進行了研究.

本文在前人研究的基礎(chǔ)上，求得了任意方向脈沖點源荷載作用下分層巖土介質(zhì)地表垂向振動位移的數(shù)值解析表達式，進行了落錘沖擊水泥公路試驗，用ANSYS/LS-DYNA進行數(shù)值模擬，使用Mathematica軟件對本文提出的方法進行編碼數(shù)值計算，將三種方法得到的地表垂向振動位移衰減曲線進行對比分析，結(jié)果可為工程中分層巖土介質(zhì)的動態(tài)響應相關(guān)計算提供一定的理論參考.

1 傳遞矩陣的推導

本文采用如圖1所示的彈性半無限空間上覆3層均勻、各向同性的分層彈性半空間巖土介質(zhì)模型，在該模型表面作用任意方向的集中脈沖點源荷載P0F(t) .

圖1 任意方向脈沖點源荷載作用下分層彈性半空間巖土介質(zhì)模型Fig.1 Layered elastic half space geotechnical model under impulsive point source load in arbitrary direction

1.1 基本方程

本文討論的情況不計體力，則分層彈性半空間巖土介質(zhì)模型任意一層內(nèi)的每一點都滿足彈性波動力學的基本方程：

1.2 積分變換

該研究問題中包含時間t 這一實變量函數(shù). 在時域中進行一些運算并不容易，但若將時間t 作拉普拉斯變換，即可將時域中的問題轉(zhuǎn)換為復數(shù)域的問題，并在復數(shù)域中作各種運算，再將運算結(jié)果作拉普拉斯反變換來求得時域中的相應結(jié)果，這就是引入拉普拉斯變換的一個主要優(yōu)點.

對公式（1）進行關(guān)于時間t的Laplace變換，得到：

對公式（4）再進行Hankel變換，可得：同理，對本構(gòu)關(guān)系式（2）和（3）也進行Laplace-Hankel聯(lián)合積分變換，可得應力的積分變換形式：

1.3 Cayley-Hamilton定理

將經(jīng)過積分變換后的基本方程和本構(gòu)方程合并，寫成常微分方程組的形式：

其中：

根據(jù)文獻［9］中關(guān)于常微分方程組的解答并結(jié)合本模型，可得i=4，故有4個方程解的形式：

可以看出，指數(shù)矩陣即描述任意深度z的狀態(tài)向量與該層表面處的狀態(tài)向量間關(guān)系的傳遞矩陣. 只要求解出該指數(shù)矩陣，就可以求得Ui(ξ,z,s) .

關(guān)于該指數(shù)矩陣的求解，根據(jù)凱利-哈密爾頓（Cayley-Hamilton）定理［10］，可以得到e[zDi(ξ,s)]的特征方程為：

求解上述方程組，得：

代入可得到傳遞矩陣：

式中：

為了方便起見，忽略式中的i，各符號意義同前. 根據(jù)推導出的傳遞矩陣，Ui(ξ,z,s)就可以寫成如下形式：

1.4 邊界條件

分層巖土介質(zhì)的層間邊界處滿足連續(xù)介質(zhì)力學的假設，所以邊界處應力位移是連續(xù)的，即：Ui(ξ,hi-1,s)=Ui-1(ξ,hi-1,s). 利用該邊界條件，再借助Ui(ξ,z,s)的公式，可建立起地表上狀態(tài)向量與分層巖土介質(zhì)及其下半無限介質(zhì)內(nèi)狀態(tài)向量的關(guān)系式：

1.5 矩陣分塊法

將邊界條件代入關(guān)系式，利用矩陣分塊法［12］求得：

化簡，得：

1.6 脈沖荷載及幅值確定

在該研究模型中，荷載的形式是任意方向的脈沖點源荷載，其形式如下：

對上式進行Laplace和0階Hankel變換：

在上述表達式中，P0為脈沖荷載的幅值. 解答在實際運用中，需確定脈沖幅值P0的大小. 通過量綱分析可知P0的單位為N. 脈沖荷載幅值考慮采用脈沖觸地瞬間的最大沖擊力. 分析對比若干種沖擊力的計算方式后，決定采用Hertz碰撞理論［13］進行求解. 觸地沖擊這一過程可看作兩彈性球體間的彈性碰撞，碰撞時的最大沖擊力P0可通過下式計算：

式中：m、M分別為兩彈性球A、B的質(zhì)量；V為碰撞前的速度；r、R分別為兩彈性球的半徑；E1、E2分別為兩彈性球的彈性模量；v1、v2為兩彈性球的泊松比. 假定沖擊體為A，地面為B，地面B的質(zhì)量M和半徑R可看作無窮大，上式可簡化為：

式（31）分母中，泊松比與彈性模量的比值大小相比較兩彈性模量的比值可忽略不計，因此最大沖擊力P0的計算公式可簡化為：

2 落錘沖擊水泥公路試驗

以落錘沖擊地面作為振源，水泥公路作為試驗對象. 為盡可能減少外界干擾對試驗的影響，同時確保試驗人員的人身安全，最終選擇南京市鼓樓區(qū)虎踞路某圍閉施工的道路區(qū)間進行. 水泥路面自上而下總共由四部分組成，依次為10 cm 厚水泥路面、10 cm 水穩(wěn)碎石基層、10 cm 三合土底基層及天然紅土基層，如圖2所示，落錘及各土層材料物理力學參數(shù)如表1所示.

圖2 水泥路面結(jié)構(gòu)Fig.2 Cement pavement structure

表1 材料物理力學參數(shù)Tab.1 Physical and mechanical parameters

2.1 試驗設備及方法

試驗的設備主要由落錘沖擊加載裝置及振動監(jiān)測裝置兩部分組成. 落錘沖擊加載裝置由三腳支撐架起吊裝置和落錘組成，振動監(jiān)測裝置采用加拿大Instantel公司生產(chǎn)的Minimate-Pro4振動和過壓監(jiān)測儀. 試驗共計取監(jiān)測點12個，監(jiān)測點編號按距沖擊點距離的遠近依次為P1～P12，每個監(jiān)測點之間的距離為0.3 m. 其中，監(jiān)測點P1距落錘沖擊點0.3 m. 按測點布置圖3將三向檢波器分別放置在對應位置，并將三向檢波器徑向指向振源方向，利用石膏將三向檢波器與地面剛性連接，待石膏凝結(jié)硬化后開始試驗，如圖4所示.

圖3 測點布置示意圖Fig.3 Layout of measuring points

圖4 三向檢波器與地面剛性連接Fig.4 The three-way geophone rigidly connected with the ground

采樣頻率設置應遵循奈奎斯特定理，即為保證離散數(shù)字信號能夠不失真地恢復到連續(xù)的模擬信號，采樣頻率應不小于信號中最高頻率的2倍. 在實際應用中，采樣頻率通常設定為信號最高頻率的5～10倍，為了更加準確地研究振動波形，采樣頻率設為信號最高頻率的10倍以上. 故試驗振動波形監(jiān)測裝置的采樣頻率設為2048 Hz. 根據(jù)現(xiàn)有試驗條件，在保證安全的前提下，落錘沖擊試驗的最終下落高度定為1.2 m，得到的振速時程曲線如圖5所示（由于測點多，此處取測點P7的波形圖用以展示）.

圖5 1.2 m沖擊高度振速時程圖（測點：P7）Fig.5 Time history diagram of vibration velocity at 1.2 m impact height（measuring point：P7）

2.2 數(shù)據(jù)處理

為了用該計算方法來驗證該試驗并與試驗結(jié)果相比較分析，我們需要得到地表垂向振動位移與距振源距離之間關(guān)系的變換曲線，所以需要對測得的速度時程曲線進行積分處理. 若直接積分，曲線會出現(xiàn)較為嚴重的漂移現(xiàn)象. 因此，還需要通過基線校正來解決該問題. 采用DPLOT 曲線處理軟件對速度時程曲線圖進行處理，通過“Generate->Integrate a Curve”指令對速度時程曲線積分得到位移時程曲線，再通過指令“Edit->Remove Trend”并選擇剛得到的位移時程曲線來對位移進行一次線性基線校正，即可得到不漂移的位移時程曲線. 利用該方法對P1～P12 監(jiān)測點的振動速度時程曲線進行處理，然后提取沿各測點位移時程曲線中的位移峰值，可得到位移隨距振源距離變換的曲線，如圖6所示.

圖6 地表垂向振動位移與振源距離關(guān)系曲線Fig.6 Relationship curve between surface vertical vibration displacement and vibration source distance

3 ANSYS/LS-DYNA數(shù)值模擬

為了與試驗結(jié)果相比對，采用ANSYS/LS-DYNA 軟件建立了4 m×1 m、厚度為10 cm 的分層巖土介質(zhì)落錘沖擊振動模型，并模擬其受2.5 kg鋼制落錘沖擊的工況. 設置總計算時間為3.6×10-3s，動畫的輸出間隔為5×10-5s，因此共輸出約75 幀. 為了讓觀測點（圖7 中黃色點所示）輸出更多的數(shù)據(jù)，采用DATABASE_HISTORY_NODE_SET 和DATABASE_NODOUT，設置觀測點的輸出間隔為1×10-6s，因此每個點可以輸出3600 個數(shù)據(jù). 在落錘的沖擊作用下，模型應力響應如圖8所示.

圖7 分層巖土介質(zhì)模型的數(shù)據(jù)觀測點Fig.7 Observation points of layered rock and soil media model

圖8 不同時刻的位移云圖Fig.8 Displacement maps at different times

根據(jù)模擬結(jié)果，選取距沖擊點0.5 m至3.8 m之間的12個監(jiān)測點P1～P12，各測點距離及垂向振動位移峰值整理如表2所示.

表2 垂向振動位移峰值Tab.2 Peak value of vertical vibration displacement

4 數(shù)值實施及計算

針對試驗的工況，利用表達式（25）進行數(shù)值計算，用來考核該解答數(shù)值化實施的可行性. 從表達式（25）可以看出：無窮域的Hankel 變換和Laplace 變換將直接影響解答數(shù)值實施的精度，所以找到相應的快速有效的計算方法將對數(shù)值化實施的可行性至關(guān)重要. 這里，選取兩種解決無窮域積分的數(shù)值計算方法，即快速Hankel 變換［14］和改進后的Durbin 方法［15-16］. 利用這兩種方法分別對表達式（25）中的無窮域Hankel變換和Laplace變換進行處理，這樣對表達式（25）進行的數(shù)值計算的實施將會變得更加方便，而且公式中矩陣的階數(shù)為4 階，因此對工作量的增加將會比較有限. 利用Mathematica 軟件編寫代碼數(shù)值計算，對結(jié)果進行數(shù)據(jù)處理，得到地表垂向振動位移與振源距離關(guān)系曲線，結(jié)合現(xiàn)場實驗的結(jié)果繪制得到對比圖（圖9）.

從圖9可以看出，利用該數(shù)值方法的計算結(jié)果與數(shù)值模擬結(jié)果更相近，二者與試驗實測值的誤差均未超過1.7 mm/s，得到的衰減規(guī)律均與試驗實測的結(jié)果衰減規(guī)律基本一致，故可證明該數(shù)值分析方法實施的可行性，結(jié)果可為相關(guān)工程計算與分析提供參考.

試驗得到的垂向振動位移衰減曲線存在明顯的起伏現(xiàn)象. 這是因為落錘作為試驗中脈沖載荷的輸入形式，實質(zhì)上可視為均布的一系列脈沖點源的集合. 一方面，對地面某測點而言，各脈沖點源相應的相位是互異的，即不同點發(fā)出的波動有相位差，振動波傳至該點合成時根據(jù)相位相同或不同分別產(chǎn)生疊加和抵消作用，致使垂向振動位移衰減曲線出現(xiàn)起伏；另一方面，應力波在層狀地基中傳播時，在交界層面將發(fā)生反射、折射，而反射波和由振源傳出的振動波發(fā)生干涉也會導致衰減曲線出現(xiàn)起伏.

圖9 地表垂向振動位移衰減曲線對比Fig.9 Comparison of attenuation curves of surface vertical vibration displacement

5 結(jié)論

采用理論模型建立、試驗研究、數(shù)值模擬及數(shù)值算法編制與程序開發(fā)相結(jié)合的研究思路，對分層巖土介質(zhì)的振動響應研究進行了系統(tǒng)深入的研究，主要成果如下：

1）在前人研究的基礎(chǔ)上，利用Laplace-Hankel 聯(lián)合積分變換和傳遞矩陣法給出了與水平面呈任意夾角的點源脈沖荷載作用下彈性半空間分層巖土介質(zhì)地表垂向位移的數(shù)值解析式.

2）進行落錘沖擊水泥公路試驗和ANSYS/LS-DYNA 軟件進行數(shù)值模擬，得到距沖擊源不同位置處地表各點的垂向振動位移，并繪制出垂向振動位移衰減曲線，與提出的數(shù)值解析方法得到的計算結(jié)果進行對比分析. 結(jié)果表明，將傳遞矩陣法應用到分層巖土介質(zhì)地表振動響應問題研究中，能較好地反映任意方向脈沖激勵下地面振動的衰減過程，且計算方法簡潔，易于應用.

3）只給出了任意方向點源脈沖荷載作用下彈性半空間分層巖土介質(zhì)地表的垂向位移表達式，對于其他荷載類型作用下的振動位移求解可以通過荷載函數(shù)與脈沖響應的卷積得到，該研究成果可為分層巖土介質(zhì)在沖擊荷載下的地表振動規(guī)律分析提供一定的指導意義.