層狀場地的精確動力剛度分析

抗波基本原理通過地震害調(diào)查和強(qiáng)震觀測，揭示了表土層附近地震特征對地震波的振幅和頻率有顯著影響。Thomson和Haskell開創(chuàng)性地給出了層狀場地中波傳播問題的傳遞矩陣方法。Kausel和Roesset利用Thomson-Haskell傳遞矩陣方法給出了層狀場地的剛度矩陣,Biot和Wolf也分別采用不同方法給出了層狀場地的剛度矩陣。尤其是,Wolf基于層狀場地動力剛度矩陣,建立了層狀場地中波的傳播以及土-結(jié)構(gòu)相互作用等問題比較完整的一套理論,并在工程中得到廣泛應(yīng)用。但上述關(guān)于剛度矩陣的研究,仍局限于二維層狀場地。本文在Wolf理論基礎(chǔ)上,給出了三維層狀場地的精確動力剛度矩陣,并利用三維層狀場地的精確動力剛度矩陣,計算分析了三維層狀場地的動力響應(yīng),給出了三維層狀場地表面或埋置矩形均布荷載及集中荷載的動力格林函數(shù)。1asv的計算,cz2mx+fmtasvix+mszcszb,bv.三維土層中的運(yùn)動如圖1所示,它是分析層狀場地的基本單元。局部坐標(biāo)系的原點(diǎn)位于厚度為d的土層頂面,z軸正向垂直向下。土層介質(zhì)為彈性均勻各向同性材料,且任一土層的材料特性為常數(shù)。以位移形式表示的穩(wěn)態(tài)三維動力平衡方程:(λ+G)(u,xx+v,xy+w,xz)+G(u,xx+v,yy+w,zz)=-ρω2u(1a)(λ+G)(u,xy+v,yy+w,yz)+G(u,xx+v,yy+w,zz)=-ρω2v(1b)(λ+G)(u,xz+v,yz+w,zz)+G(u,xx+v,yy+w,zz)=-ρω2w(1c)設(shè)滿足方程(1)的P波和S波的波勢函數(shù)分別為:e=-iωcpApexp[iωcp(-lxx-lyy-lzz)](l2x+l2y+l2z=1)(2)Ω=-iω2csCexp[iωcs(-mxx-myy-mzz)](m2x+m2y+m2z=1;mxCx+myCy+mzCz=0)(3)由{u,v,w}Τ=?e+?×Ω,同時進(jìn)一步將S波的位移分量分解成為幅值為ASH的水平分量,以及在由整體坐標(biāo)軸z與波的傳播方向所決定的平面內(nèi)幅值為ASV的另一分量,可求得位移:u(kx,ky,z)=lxApexp(-iω(lxcpx+lycpy+lzcpz))+mxmz√m2x+m2yASVexp(-iω(mxcsx+mycsy+mzcsz))-my√m2x+m2yASΗexp(-iω(mxcsx+mycsy+mzcsz))(4a)v(kx,ky,z)=lyApexp(-iω(lxcpx+lycpy+lzcpz))+mymz√m2x+m2yASVexp(-iω(mxcsx+mycsy+mzcsz))+mx√m2x+m2yASΗexp(-iω(mxcsx+mycsy+mzcsz))(4b)w(kx,ky,z)=lzApexp(-iω(lxcpx+lycpy+lzcpz))-√m2x+m2yASVexp(-iω(mxcsx+mycsy+mzcsz))(4c)其中,ASΗ=Cz√m2x+m2y,ASV=(mxCy-myCx)√m2x+m2y。考慮到必須滿足所有的6個邊界條件(幅值為u、v和w的位移,幅值為τzx、τzy和σz的應(yīng)力),土層中必須包含傳播方向與z軸正向相反的來波和傳播方向與z軸正向相同的去波,引入材料阻尼后,土層中的位移可表示為:u(kx,ky,z)=lx[Apexp(iωlzc*pz)+Bpexp(-iωlzc*pz)]exp(-iω(lxc*px+lyc*py))-mxmz√m2x+m2y[ASVexp(iωmzc*sz)-BSVexp(-iωmzc*sz)]exp(-iω(mxc*sx+myc*sy))-my√m2x+m2y[ASΗexp(iωmzc*sz)+BSΗexp(-iωmzc*sz)]exp(-iω(mxc*sx+myc*sy))(5a)v(kx,ky,z)=ly[Apexp(iωlzc*pz)+Bpexp(-iωlzc*pz)]exp(-iω(lxc*px+lyc*py))-mymz√m2x+m2y[ASVexp(iωmzc*sz)-BSVexp(-iωmzc*sz)]exp(-iω(mxc*sx+myc*sy))+mx√m2x+m2y[ASΗexp(iωmzc*sz)+BSΗexp(-iωmzc*sz)]exp(-iω(mxc*sx+myc*sy))(5b)w(kx,ky,z)=-lz[Apexp(iωlzc*pz)-Bpexp(-iωlzc*pz)]exp(-iω(lxc*px+lyc*py))-√m2x+m2y[ASVexp(iωmzc*sz)+BSVexp(-iωmzc*sz)]exp(-iω(mxc*sx+myc*sy))(5c)其中,lx、ly、lz為P波傳播方向的方向余弦,lx=cos(φαp)、ly=cos(φβp)、lz=cos(φγp),而φαp、φβp和φγp為P波傳播方向的方向角;mx、my、mz為S波傳播方向的方向余弦,mx=cos(ψαs)、my=cos(ψβs)、mz=cos(ψγs),而ψαs、ψβs和ψγs為S波傳播方向的方向角;Ap、ASV和ASH分別為來波P、SV和SH波的幅值;Bp、BSV和BSH分別為去波P、SV和SH波的幅值;λ*=λ(1+2iζ),G*=G(1+2iζ),ζ為材料滯回阻尼比;c*p=√(λ*+2G*)ρ,c*s=√G*ρ,ρ為材料質(zhì)量密度。此種形式的位移使得土層的頂面和底面的邊界條件也按exp(-iωlxx/c*p)和exp(-iωlyy/c*p)與exp(-iωmxx/c*s)和exp(-iωmyy/c*s)形式發(fā)生變化。為了使它們分別隨x和y有相同的變化規(guī)律,條件lx/c*p=mx/c*s和ly/c*p=my/c*s必須滿足。記x方向的波數(shù)kx=lxω/c*p=mxω/c*s,y方向的波數(shù)ky=lyω/c*p=myω/c*s,z方向的縱波和橫波波數(shù)分別為kpz=lzω/c*p,ksz=mzω/c*s,這樣,式(5)可寫為:u=u(kx,ky,z)exp(-ikxx)exp(-ikyy)(6a)v=v(kx,ky,z)exp(-ikxx)exp(-ikyy)(6b)w=w(kx,ky,z)exp(-ikxx)exp(-ikyy)(6c)由τzx=G*(u,z+w,x),τzy=G*(v,z+w,y)和σz=λ*(u,x+v,y+w,z)+2G*w,z,可得:τzx=C11Apexp(ikpzz)+C12Bpexp(-ikpzz)+C13ASVexp(ikszz)+C14BSVexp(-ikszz)+C15ASΗexp(ikszz)+C16BSΗexp(-ikszz)(7a)τzy=C21Apexp(ikpzz)+C22Bpexp(-ikpzz)+C23ASVexp(ikszz)+C24BSVexp(-ikszz)+C25ASΗexp(ikszz)+C26BSΗexp(-ikszz)(7b)σz=C31Apexp(ikpzz)+C32Bpexp(-ikpzz)+C33ASVexp(ikszz)+C34BSVexp(-ikszz)+C35ASΗexp(ikszz)+C36BSΗexp(-ikszz)(7c)由于在z為常數(shù)的截面上無另外三個應(yīng)力(τxy、σx和σy)作用,故對此不予計算。由式(6)和式(7),分別令z=0和d,可得在局部坐標(biāo)系下土層頂面位移u1、v1、w1和應(yīng)力τzx1、τzy1、σz1,以及土層底面的位移u2、v2、w2和應(yīng)力τzx2、τzy2、σz2;它們均分別為土層中來波和去波的幅值A(chǔ)p、Bp、ASV、BSV、ASH、BSH的函數(shù),為使矩陣對稱,將z向的位移和應(yīng)力乘以i:[u1v1iw1u2v2iw2]=[lxlx-mxmzβ5mxmzβ5-myβ5-myβ5lyly-mymzβ5mymzβ5mxβ5mxβ5-ilzilz-i/β5-i/β500lxβ1lxβ2-mxmzβ3β5mxmzβ4β5-myβ3β5-myβ4β5lyβ1lyβ2-mymzβ3β5mymzβ4β5mxβ3β5mxβ4β5-ilzβ1ilzβ2-iβ3/β5-iβ4/β500][ApBpASVBSVASΗBSΗ](8a)[-τzx1-τzy1-iσz1τzx2τzy2iσz2]=[C11C12C13C14C15C16C21C22C23C24C25C26C31C32C33C34C35C36C11β1C12β2C13β3C14β4C15β3C16β4C21β1C22β2C23β3C24β4C25β3C26β4C31β1C32β2C33β3C34β4C35β3C36β4][ApBpASVBSVASΗBSΗ](8b)其中:β1=exp(ikpzz);β2=exp(-ikpzz);β3=exp(ikszz);β4=exp(-ikszz);β5=1√m2x+m2y;C11=-i(lxkpz+lzkx)G*;C12=i(lxkpz+lzkx)G*;C13=i(mxmzksz-(m2x+m2y)kx)G*β5;C14=i(mxmzksz-(m2x+m2y)kx)G*β5;C15=imykszG*β5;C16=-imykszG*β5;C21=-i(lykpz+lzky)G*;C22=i(lykpz+lzky)G*;C23=i(mymzksz-(m2x+m2y)ky)G*β5;C24=i(mymzksz-(m2x+m2y)ky)G*β5;C25=-imxkszG*β5;C26=imxkszG*β5;C31=-(lxkx+lyky+lzkpz)λ*-2lzkpzG*;C32=-(lxkx+lyky+lzkpz)λ*-2lzkpzG*;C33=(λ*(mxmzkx+mymzkx)-(λ*+2G*)(m2x+m2y)ksz)β5;C34=((λ*+2G*)(m2x+m2y)ksz-λ*(mxmzkx+mymzkx))β5;C35=(mykx-mxky)λ*β5;C36=(mykx-mxky)λ*β5。將這些波的幅值消去,可得土層的動力剛度矩陣。在單元的負(fù)微元面(頂面)上,應(yīng)力τzx1、τzy1、σz1與坐標(biāo)系方向相反為負(fù)。在單元的正微元面(底面)上,應(yīng)力τzx2、τzy2、σz2與坐標(biāo)系方向相同為正。則在整體坐標(biāo)系下定義的外荷載幅值和在局部坐標(biāo)系下定義的應(yīng)力幅值關(guān)系為:Px1=-τzx1;Py1=-τzy1;Pz1=-σz1;Px2=τzx2;Py2=τzy2;Pz2=σz2。最后可得土層的動力剛度矩陣為(矩陣中各元素詳見附錄):{Ρx1,Ρy1,iΡz1,Ρx2,Ρy2,iΡz2}Τ=SLΡ-SV-SΗ{u1,v1,iw1,u2,v2,iw2}Τ(9)2半空間動力剛度矩陣在半空間表面施加荷載只會產(chǎn)生幅值為Bp、BSV和BSH的去波,由于必須滿足輻射條件(波從無窮遠(yuǎn)處傳至自由表面的能量為零),半空間中不包含幅值為Ap、ASV和ASH的來波。以上標(biāo)R表示半空間,在式(6)和(7)中令A(yù)p、ASV、ASH和z為0,Px0=-τzx1;Py0=-τzy1;Pz0=-σz1,消去Bp、BSV和BSH后得半空間動力剛度矩陣SRP-SV-SH(矩陣中各元素詳見附錄):{Ρx0,Ρy0,iΡz0}Τ=SRΡ-SV-SΗ{u0,v0,iw0}Τ(10)與二維情況相同,對于無阻尼介質(zhì),SRP-SV-SH是虛數(shù),它相應(yīng)于一個阻尼器;相反,SLP-SV-SH是實(shí)數(shù),可理解成一個彈簧。對于半空間,能量向無窮遠(yuǎn)輻射。在有阻尼體系中,令d→∞,則土層動力剛度矩陣對角線上的子矩陣(3×3階)收斂于半空間的動力剛度矩陣,而非對角線矩陣趨于零。3整體動力剛度矩陣圖2表示一n個土層置于半空間上。由自由表面零應(yīng)力條件和各土層交界面的連續(xù)條件,集整土層和半空間的剛度矩陣SLP-SV-SH和SRP-SV-SH,可得整體動力剛度矩陣SΡ-SV-SΗ,則場地動力平衡方程為:SΡ-SV-SΗU=Q(11)式中,U為位移列向量,由u1、v1、iw1…uj、vj、iwj…un+1、vn+1、iwn+1組成;Q是外荷載向量,當(dāng)作用外荷載時,與外荷載相應(yīng)位置的元素不為0,其他元素則為0。計算自由場時,若控制點(diǎn)選在基巖露頭,則Q中最后三個元素由式(10)確定,其他元素為0。4自由場地表位移與文獻(xiàn)的對比對均勻半空間,由式(7)引入自由表面邊界條件(τzx=0,τzy=0,σz=0;z=0),則可得用入射波幅值A(chǔ)Rp、ARSV和ARSH來表示的反射波幅值BRp、BRSV和BRSH,然后再代入式(6)便可得地表位移u0、v0和w0與ARp、ARSV和ARSH的關(guān)系;對基巖上層狀場地,將控制點(diǎn)選在基巖露頭,記u0、v0和w0為基巖露頭運(yùn)動,則Q中最后三個元素可由式(10)確定,其他元素為0,然后利用直接剛度法,由式(11)可得整個場地的運(yùn)動。為了驗(yàn)證本文的正確性,圖3給出了當(dāng)入射波在x-z平面內(nèi)(ly=my=0)時的自由場地表位移與文獻(xiàn)的對比。三個圖分別為:(a)無阻尼均勻半空間在入射SV波時,不同泊松比情況地表位移隨入射角度的變化;(b)基巖(半空間)上單一土層,對應(yīng)不同角度入射SV波,地表位移隨入射波頻率的變化情況(計算參數(shù):基巖與土層的剪切波速比為5.0、質(zhì)量密度比為1.0,基巖和土層的阻尼比均為0.05,泊松比均為0.33);(c)基巖(半空間)上單一土層,對應(yīng)不同角度SV波,土層底面與頂面位移幅值比隨入射頻率的變化情況(計算參數(shù):基巖與土層的剪切波速比為5.0、密度比為1.0,基巖和土層的阻尼比均為0.05,泊松比均為0.33)。從圖可以看出,當(dāng)入射波在某一坐標(biāo)平面內(nèi)時,自由場位移與二維情況完全相同,說明三維土層和半空間剛度矩陣可完全退化為文獻(xiàn)二維土層和半空間剛度矩陣。圖4給出了基巖(半空間)上單一土層,分別入射P、SV和SH波時地表位移幅值隨入射頻率的變化情況。基巖與土層的剪切波速比cRs/cLs=5.0、質(zhì)量密度比ρR/ρL=1.0、基巖和土層的阻尼比分別為ζL=0.05和ζR=0.02、基巖和土層的泊松比νR=νL=0.33。從圖中可以看出,地表位移幅值隨著入射頻率的變化而變化,且出現(xiàn)多個峰值;另外,當(dāng)入射SH波時,地表只有水平位移。5b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b-b0.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.4.3.3.4.3.3.4.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.4.3.4和3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.作為剛度矩陣的另一重要應(yīng)用,下面給出表面和埋置矩形均布荷載以及集中荷載的動力格林影響函數(shù)的計算方法。如圖5所示。對于矩形均布荷載px0、py0和pz0,沿兩個水平方向展開成關(guān)于exp(-ikxx)exp(-ikyy)的傅里葉積分,得到(kx,ky)域上的荷載幅值:px(kx,ky)=px04π2∫Δa-Δa∫Δb-Δbexp(ikxx)exp(ikyy)dxdy=px0sin(kxΔa)sin(kyΔb)π2kxky(12a)py(kx,ky)=py04π2∫Δa-Δa∫Δb-Δbexp(ikxx)exp(ikyy)dxdy=py0sin(kxΔa)sin(kyΔb)π2kxky(12b)pz(kx,ky)=pz04π2∫Δa-Δa∫Δb-Δbexp(ikxx)exp(ikyy)dxdy=pz0sin(kxΔa)sin(kyΔb)π2kxky(12c)而對于表面或埋置集中荷載情況,亦即,Δa→0Δb→0,上述解答可退化為:px(kx,ky)=limΔa,Δb→0px0sin(kxΔa)sin(kyΔb)π2kxky=px0ΔaΔbπ2=Ρx04π2(13a)py(kx,ky)=limΔa,Δb→0py0sin(kxΔa)sin(kyΔb)π2kxky=py0ΔaΔbπ2=Ρy04π2(13b)pz(kx,ky)=limΔa,Δb→0pz0sin(kxΔa)sin(kyΔb)π2kxky=pz0ΔaΔbπ2=Ρz04π2(13c)(kx,ky)域上的位移幅值u可由式(11)計算。外荷載列向量Q中Q3(j-1)+1=px(kx,ky)、Q3(j-1)+2=py(kx,ky)、Q3j=ipz(kx,ky),其他元素為0,j為作用荷載的土層。由各層的位移和式(6)可得各土層來波和去波的幅值系數(shù),然后由式(6)和式(7)可得場地任意點(diǎn)的位移和應(yīng)力,最后進(jìn)行傅里葉逆變換可得空間域內(nèi)場地任意點(diǎn)的位移和應(yīng)力幅值:F(x,y,z)=∫+∞-∞∫+∞-∞F(kx,ky,z)exp(-ikxx)exp(-ikyy)dkxdky(14)其中,F(x,y,z)和F(kx,ky,z)分別為空間域和波數(shù)域的位移或應(yīng)力幅值。5.1各因素作用下的位移和激勵頻率為驗(yàn)證本文方法的精度,表1和表2分別給出了均勻半空間內(nèi)部豎向靜態(tài)集中荷載(以非常低的頻率模擬靜力問題:無量綱激勵頻率η=ωa/cs=0.001)作用下無量綱位移wGz0pz0和無量綱應(yīng)力σzzz20pz0與文獻(xiàn)和結(jié)果的對比。可以看出,無論是位移還是應(yīng)力均與文獻(xiàn)吻合很好。表3給出了簡諧正方形均布荷載(邊長為2a)作用下均勻半空間表面正方形中心的位移uG/(px0a2)和wG/(pz0a2)與文獻(xiàn)的對比,無量綱激勵頻率η=ωa/cs。可以看出,本文結(jié)果與文獻(xiàn)吻合也很好。5.2不同激勵頻率的時域演化以均勻半空間(無阻尼)和基巖上單一土層(有阻尼)為例,計算埋置正方形均布荷載(2Δa=2Δb=a,面積S=a2)產(chǎn)生的位移和應(yīng)力。土層厚度h/a=3.0,荷載作用深度為z0/a=1.0;對于無阻尼均勻半空間,阻尼比ζ=0.001,