熊 浩 高廣運 馮世進 歐湘萍

(臺州學(xué)院建筑工程學(xué)院1) 臺州 318000) (同濟大學(xué)地下建筑與工程系2) 上海 200092)

(同濟大學(xué)巖土及地下工程教育部重點實驗室3) 上海 200092) (武漢理工大學(xué)交通學(xué)院4) 武漢 430063)

0 引 言

彈性波傳播的數(shù)值模擬在鐵道及城市軌道交通環(huán)境振動分析和隔振設(shè)計、地震勘探等領(lǐng)域有廣泛的應(yīng)用.有限差分法、有限元法和偽譜法是三種主要的數(shù)值模擬方法.Virieux等[1]提出了基于應(yīng)力-速度混合變量彈性波方程的交錯網(wǎng)格差分法,提高了差分法的計算精度和穩(wěn)定性.裴正林[2]將交錯網(wǎng)格計算技術(shù)與高精度差分計算格式結(jié)合起來,用于三維波動模擬.Chen[3]則進一步提出了旋轉(zhuǎn)交錯網(wǎng)格計算技術(shù),增加了差分法模擬非均勻介質(zhì)中波動問題的能力.有限元法雖理論成熟,適應(yīng)性強,但計算所需存儲空間較大,計算效率稍差.偽譜法[4]則采用快速傅立葉變換計算空間導(dǎo)數(shù),用差分法計算時間導(dǎo)數(shù).通過引入映射函數(shù)雖可以使其在計算中考慮地表形狀的變化,但如何合理選擇映射函數(shù)有一定困難.

近年來,張劍鋒[5-6]給出了基于應(yīng)力-速度混合變量彈性波方程及任意四邊形網(wǎng)格差分算子的數(shù)值方法——格子法,它融合了有限元可以適應(yīng)復(fù)雜形狀邊界條件、有限差分法計算簡便及無須計算剛度矩陣的優(yōu)點,可以方便地引入任意形狀自由表面邊界條件.

本文基于以位移為變量彈性波方程,給出了二維格子法的計算方法,推導(dǎo)適合該法的人工邊界條件,并通過經(jīng)典Lam b問題算例驗證了該方法的正確性.與傳統(tǒng)的基于混合變量的格子法相比,該方法采用單點高斯積分法計算線積分,提高了積分的準(zhǔn)確性;文中所給出的人工邊界條件計算簡便,透射外行波效果良好.應(yīng)用上述數(shù)值計算方法,本文分析了在不同填充溝深度條件下的隔振效果,即地表豎向振幅衰減系數(shù)變化.

1 計算方法

1.1 基本計算公式

在笛卡爾坐標(biāo)系中,不計體力,均質(zhì)彈性介質(zhì)中彈性波傳播控制方程為

式中:ρ為彈性介質(zhì)的密度,kg/m3;μ和λ為介質(zhì)的Lamé系數(shù),Pa;σx,σy,σz和 τxy為應(yīng)力張量,Pa;ux和uy為位移,m.

將計算區(qū)域以三角形或四邊形網(wǎng)格進行離散,局部網(wǎng)格如圖1所示.圖中的實線是離散時所生成網(wǎng)格線,實心點j,k,l,…是網(wǎng)格線頂點;圖中的虛線是各段實線中點與網(wǎng)格中心點的連線,四邊形網(wǎng)格中空心點(如1,2,3,…)是各段虛線的中點,而三角形網(wǎng)格中的空心點(如點7)是三角形的中心點.如果在實心點處定義位移分量ui(i=x,y),那么對于四邊形網(wǎng)格而言,可以利用任意四邊形網(wǎng)格差分算子[6]計算位移在域內(nèi)任意點處的空間一階導(dǎo)數(shù)值.

圖1 局部網(wǎng)格示意圖

對于三角形網(wǎng)格,以圖1中三角形 fm j為例,其中心點處的位移空間一階導(dǎo)數(shù)計算公式為[7]

式中:A為三角形 fm j的面積;(ui)f,(ui)m和(ui)j為三個頂點f,m,j處的位移值;bf=(yjym)/2,cf=(xj-xm)/2,(f→m→j→t),通過下標(biāo)輪換可以類似地計算出 bm,cm等;(xj,yj),(xm,ym)分別為j,m點處的坐標(biāo).

從圖1可見,各個實心點周圍的虛線都圍成一個閉合區(qū)域.假定結(jié)點 f周圍虛線所圍的區(qū)域為 Ω,其邊界1-2-3-4-5-6-7記為 Γ.在 Ω上對式(1)積分可得

式中:α,β分別為邊界 的外法線方向余弦.采用集中質(zhì)量模型,式(4)～(5)中的等號左端的面積積分可簡化為僅與 f點加速度有關(guān),即分別為.式中:M f是區(qū)域 Ω的質(zhì)量.對于式(3)右端的線積分,在邊界Γ每一段上采用單點高斯積分公式進行數(shù)值計算.此時,每段上的積分點正是邊界 Γ上的各個空心點.設(shè)式(4)和式(5)線積分值分別為(Ax)f和(Ay)f,則

對式(6)中二階導(dǎo)數(shù)應(yīng)用中心差分公式計算,可得

式中:上標(biāo) n為nΔt時刻,類似地n+1代表(n+1)Δt時刻,n-1為(n-1)Δt時刻;Δt為計算時間步長.根據(jù)式(7),可求得

1.2 計算步驟

對節(jié)點 f而言,應(yīng)用單點高斯積分公式,在邊界Γ上計算式(3)中線積分(Ax)f和(Ay)f,并把它們代入到式(7)中,計算出(n+1)Δt時刻的位移(un+1x)f,(un+1y)f.這樣,就完成了 f點從nΔt時刻到(n+1)Δt時刻的位移分量遞推計算.對于其他實心點,可以類似于 f點進行處理.從而整個實線網(wǎng)格上的全部節(jié)點在(n+1)Δt時刻位移可以完全確定.類似地,可以計算出(n+1)Δt時刻的位移場,直到設(shè)定的計算時間終點為止.

為了保證計算的穩(wěn)定,時間步長Δt與空間最小網(wǎng)格邊長h應(yīng)滿足穩(wěn)定條件:Δt≤h/vp,式中:v p為縱波的波速.

1.3 自由表面邊界條件

假定在圖1中h,l,k是自由表面邊界上的節(jié)點.在 l點處,在圍線 l-Ma-9-3-2-8-Mb所圍區(qū)域內(nèi)類似于式(4)～(5)建立平衡方程

式中:M l為上述圍線所圍區(qū)域的質(zhì)量.

式(9)中線積分分為兩部分,其中Ma-9-3-2-8-Mb線段上的積分計算方法與前述式(4)～(5)中的線積分計算辦法完全相同.而Mbl-Ma線段上的積分事實上就等于作用于表面邊界上的已知外荷載.因此,在自由邊界上邊界條件可以方便地引入到計算中.

2 人工邊界條件

考慮到實際上波的傳播是在半無限域中進行的,而計算模型只能是有限的,因此必須在計算區(qū)域的邊界上引入人工邊界條件,以準(zhǔn)確地模擬波的無限傳播.本文給出一種簡單實用的人工邊界條件.

將式(2)代入到式(1)中,可得

假定xoy坐標(biāo)系旋轉(zhuǎn)θ角得到新坐標(biāo)系xoy′,其坐標(biāo)軸x′與平面上某處波的傳播方向相一致,忽略另一個軸方向波的傳播,亦即假定?/?y′=0.兩個坐標(biāo)系之間有如下關(guān)系式

由此,可得

把式(13)代入到式(10～11)中,可得以矩陣表示的方程

對Q作特征值分解,Q=P-1ΛP,于是式(14)可以寫成

把xoy系中PU表達式代入上式中,可得

式中 :ux(vs)=ux(x-vsΔt cosθ,y+vsΔt sin θ,t-Δt),其余類似.由式(18)不難得到ux,yy的計算表達,可以發(fā)現(xiàn)邊界節(jié)點位移由前一時刻的鄰近節(jié)點的位移計算確定.在不同的邊界上,應(yīng)選擇不同的θ,例如在邊界法線與 x正方向平行的邊界上,應(yīng)取θ=0°計算,而在邊界法線方向與y軸的負方向上一致的邊界上,則應(yīng)取θ=90°計算.

3 計算方法的驗證

為了本文計算方法的正確性,這里計算了Lamb問題,并將計算結(jié)果與解析解[8]進行了對比.假定半無限地基密度 ρ=2 250 kg/m3,其Lamé系數(shù) μ=1.048×103MPa,λ=1.863×103MPa;計算區(qū)域為2 000 m×1 000 m的矩陣區(qū)域,采用邊長為5 m的正方形網(wǎng)格,時間步長取為0.002 s;振源為一垂直作用在地表中點的集中荷載,其表達式為 f(t)=exp(-200(t-0.178)2);計算采用上述人工邊界條件.圖2同時給出了地表上距垂直力源距離為500 m的點水平位移的解析解和本文數(shù)值計算結(jié)果,并給出了無人工邊界時的數(shù)值解結(jié)果.由圖可知,采用本文所述的計算方法及人工邊界條件所得的數(shù)值結(jié)果與解析解答基本吻合,由此可見本文方法的正確性;對比無人工邊界條件計算結(jié)果與有人工邊界條件計算結(jié)果,可見本文所述的人工邊界可以較好地透射出外行波,是一種有效的人工邊界條件.

圖2 水平位移分量的解析解與數(shù)值解的比較

4 填充溝隔振分析

振動對人們的生活與工作環(huán)境造成了不可忽視的影響,振動已被列七大公害之一.治理環(huán)境振動的有效途徑之一是設(shè)置屏障隔振,填充溝就是屏障隔振措施之一.利用上述的計算方法,本文研究了采用填充溝進行隔振時,不同的填充溝深度條件下的豎向振幅衰減變化.

圖3是計算模型.均勻彈性半空間地基的密度為ρs=2 000 kg/m3;剪切模量為G s=1.923×107Pa;泊松比v s=0.3.填充溝內(nèi)填充物為橡膠,其密度ρ=1 200 kg/m3;剪切模量為G=9.655×105Pa;泊松比υ=0.48.振源是頻率為5 H z的正弦力,表達式為 f(t)=sin(10πt).假定填充溝的寬度不變,為2m;深度h變化,分別取10,12,15,18和20m;填充溝距振源為L1=50m;在填充溝后的地表,選取與填充溝距離為L2的若干個觀測點.

圖3 填充溝計算示意圖

Woods[9]研究屏障隔振時,提出在頻域分析中,對于隔振效果可采用振幅衰減系數(shù) AR來衡量.類似地,本文定義豎向振幅衰減系數(shù)

用于評價隔振效果.根據(jù)這一定義,如果ˉAR值越小,則表明隔振的效果越好,反之則表明隔振效果不佳.

圖4是計算結(jié)果.為便于分析,圖中長度采用量綱一的量的長度:L*2=L2/λR,h*=h/λR,λR 是R 波波長,λR=18.188 m.由圖可知,對于所設(shè)置的不同填充溝深度,都可以取得一定的隔振效果.在計算模型中,隨著溝深度的逐步增加,各點隔振效果先有不同程度地降低,但當(dāng)h*≈λR時,各點隔振效果卻達到最佳水平.這表明,增加填充溝深度最終是有助于增加隔振效果的.對于某一填充溝深度而言,填充溝后地表各觀測點的豎向振幅衰減系數(shù)均在 L*2=1.10λR～3.10λR,L*2=3.50 λR～4.75λR范圍內(nèi)相對偏大,而在 L*2=3.30λR和L*2=4.95λR處相對較小.這說明:一方面,填充溝的深度達一定值后其變化對于豎向振幅衰減系數(shù)ˉA R影響不大;另一方面,填充溝的隔振效果是隨距離L*2而不斷變化的,在不同的距離處隔振效果相差較大,且存在著某個點,對應(yīng)著最小豎向振幅衰減系數(shù).所以在設(shè)計填充溝時,應(yīng)注意選擇填充溝的位置,確保隔振保護區(qū)域為最小豎向振幅衰減系數(shù),以達到最佳的隔振效果.

圖4 不同填充溝深度h*時豎向振幅衰減系數(shù)隨距離L*2變化曲線

5 結(jié) 論

1)本文基于以位移為變量彈性波方程,給出了二維格子法的計算方法,并推導(dǎo)適合該法的人工邊界條件.

2)經(jīng)典的Lam b算例表明了本文計算方法的正確性和邊界條件的適應(yīng)性.該方法能很好地適應(yīng)不規(guī)則的邊界條件,且計算簡便、準(zhǔn)確,可用于彈性波傳播及散射問題研究.

3)填充溝屏障隔振分析表明,增加填充溝深度有助于增加隔振效果,但填充溝深度并不影響最小豎向振幅衰減系數(shù)點的位置.

[1]V irieux J.P-SV w ave p ropagation in heterogeneous media:velocity-stress finite-difference method[J].Geophysics,1986,52(4):888-901.

[2]裴正林.三維各向同性介質(zhì)彈性波方程交錯網(wǎng)格高階有限差分法模擬[J].石油物探,2005,44(4):308-316.

[3]Chen H,W ang X M,Zhao H B.A rotated staggered grid finite-difference with the absorbing boundary condition of a perfec tly matched layer[J].Chinese Science Bulletin,2006,51(19):2304-2314.

[4]Tessmer E,Kosloff D,Behle A.Elastic wave propagation simu lation in the p resence of surface topography[J].Geophysical Journal International.1992,108:621-632.

[5]Zhang JF,Liu T.P-SV-wave propagation in heterogeneousmedia:G ridmethod[J].Geophysical Journal International.1999,136,431-438.

[6]張劍鋒.彈性波數(shù)值模擬的非規(guī)則網(wǎng)格差分法[J].地球物理學(xué)報,1998,41(增刊):357-366.

[7]Cook R D.Concepts and applications of finite element analysis[M].New York:John Wiley&Sons,1974.

[8] Lamb H.On the propagation of tremors over the surface of an elastic so lid[J].Philosophical Transactions o f the Royal Society of London,1904,203:1-44.

[9]Woods R D,Barnet N E,Sangesser R H.A new tool for soil dynamics[J].Journal of Geotechnical Engineering,1974,100(11):1231-1247.