基于隨機(jī)激勵(lì)的土非線性隨機(jī)場(chǎng)數(shù)值模擬方法

基于隨機(jī)激勵(lì)的土非線性隨機(jī)場(chǎng)數(shù)值模擬方法

0基于理論分析思路的地震動(dòng)空間變化規(guī)律研究方法地震空間變化規(guī)律對(duì)大范圍內(nèi)的地震響應(yīng)有重要影響?；诿芗卣鹋_(tái)陣記錄的統(tǒng)計(jì)回歸工作極大地促進(jìn)了地震動(dòng)空間變化規(guī)律的研究,但是這種統(tǒng)計(jì)方法所得到的模型強(qiáng)烈依賴于密集臺(tái)陣記錄,無(wú)法有效地將其研究結(jié)果應(yīng)用于其他實(shí)際工程場(chǎng)址。為了改變這種模型研究與工程應(yīng)用之間的脫節(jié)現(xiàn)狀,研究者們開(kāi)始尋求基于理論分析思路的地震動(dòng)空間變化規(guī)律研究途徑。例如:文獻(xiàn)利用數(shù)值方法研究了從基巖垂直輸入隨機(jī)激勵(lì)下隨機(jī)介質(zhì)場(chǎng)地的地震動(dòng)相干函數(shù);文獻(xiàn)、利用波動(dòng)有限元分析技術(shù)分析了復(fù)雜場(chǎng)地地震動(dòng)空間相關(guān)性。上述研究都是在線彈性范圍內(nèi)進(jìn)行的,然而,在地震動(dòng)過(guò)程中,土介質(zhì)中的剪切波速往往出現(xiàn)顯著的降低現(xiàn)象,體現(xiàn)出強(qiáng)烈的非線性性質(zhì),因此,有必要尋求一種能夠考慮場(chǎng)地介質(zhì)非線性特性的場(chǎng)地地震動(dòng)相干函數(shù)的分析方法。1分析1.1考慮恢復(fù)力的隨機(jī)非線性不帶乘子法分析由于地震動(dòng)具有隨機(jī)性,應(yīng)采用隨機(jī)振動(dòng)的分析方法,結(jié)構(gòu)非線性隨機(jī)振動(dòng)的基本微分方程一般具有如下形式:g(¨x,˙x,x)=f(t)(1)g(x¨,x˙,x)=f(t)(1)式中x為位移向量;g(¨xx¨,˙xx˙,x)中第i分量gi(¨xx¨,˙xx˙,x)為與第i自由度相應(yīng)的非線性恢復(fù)力關(guān)系;f(t)為隨機(jī)輸入。在非線性隨機(jī)振動(dòng)分析中,目前使用較廣的是隨機(jī)等效線性化方法,該方法假定存在一個(gè)與式(1)等效的線性系統(tǒng):Μ¨x+C˙x+Κx=f(t)(2)Mx¨+Cx˙+Kx=f(t)(2)系統(tǒng)中各參數(shù)的取值應(yīng)遵循使等效線性系統(tǒng)的反應(yīng)x盡可能接近原非線性系統(tǒng)的反應(yīng)的原則?；谏鲜鏊枷?研究者一般采用嚴(yán)格的隨機(jī)等效線性化方法和基于確定性等效線性化分析方法提出的隨機(jī)等效線性化方法對(duì)問(wèn)題進(jìn)行求解。其中,嚴(yán)格的隨機(jī)等效線性化方法分析繁雜,難以推廣,在土動(dòng)力分析中較少應(yīng)用。相對(duì)而言,基于確定性等效線性化分析方法提出的隨機(jī)等效線性化方法雖然理論上并不嚴(yán)密,且其中參數(shù)如γmax的確定也帶有一定的經(jīng)驗(yàn)性,但該方法簡(jiǎn)單易行、適用范圍廣。本文把工程波動(dòng)理論中散射問(wèn)題的處理方法與隨機(jī)振動(dòng)的虛擬激勵(lì)方法相結(jié)合,提出一種考慮土介質(zhì)非線性性質(zhì)的隨機(jī)波動(dòng)分析方法,并將其應(yīng)用于一致激勵(lì)下工程場(chǎng)地地震動(dòng)相干性的分析。1.2人工邊界的透射以圖1的直入射情況為例說(shuō)明本文的計(jì)算方法。設(shè)二維波動(dòng)分析問(wèn)題,入射波為SV波,對(duì)工程場(chǎng)地進(jìn)行波動(dòng)有限元分析時(shí),從無(wú)限域中截取一個(gè)有限計(jì)算區(qū)域,在計(jì)算區(qū)域的邊界處需要采用特殊的技術(shù)模擬外行波場(chǎng),這種經(jīng)人為引入的邊界一般稱之為人工邊界。對(duì)人工邊界范圍內(nèi)的區(qū)域進(jìn)行有限元離散后,可以寫(xiě)出運(yùn)動(dòng)控制方程的波動(dòng)有限元矩陣形式:Μ¨u(t)+C˙u(t)+Cu(t)=R(t)(3)Mu¨(t)+Cu˙(t)+Cu(t)=R(t)(3)其中M,C,K,R(t)分別是單元的質(zhì)量矩陣、阻尼矩陣、剛度矩陣和荷載矢量。采用文獻(xiàn)發(fā)展的多次透射人工邊界,有如下以人工邊界范圍內(nèi)部結(jié)點(diǎn)位移表示的人工邊界節(jié)點(diǎn)位移關(guān)系式:up+10=Ν∑j=1(-1)j+1CΝj(up+1-jj-up+1-jr,j)+up+1r,0(4)up+10=∑j=1N(?1)j+1CNj(up+1?jj?up+1?jr,j)+up+1r,0(4)式中N為人工邊界的透射階數(shù),本文中N=2;CΝjNj為二項(xiàng)式組合系數(shù);up+1-jjp+1?jj表示與所考慮人工邊界結(jié)點(diǎn)處于同一法線上、且距離為jcaΔt的內(nèi)部質(zhì)點(diǎn)在(p+1-j)Δt時(shí)刻的波動(dòng)描述量;Δt為計(jì)算時(shí)步,ca為人工波速,本文取介質(zhì)的剪切波速;upr,jpr,j表示與所考慮人工邊界結(jié)點(diǎn)處于同一法線上、且距離為jcaΔt的內(nèi)部質(zhì)點(diǎn)在pΔt時(shí)刻的射入人工邊界的波的波動(dòng)描述量,正確確定由射入人工邊界的波在人工邊界范圍內(nèi)形成的波場(chǎng)是波動(dòng)有限元的關(guān)鍵問(wèn)題,入射波場(chǎng)的計(jì)算方法見(jiàn)文獻(xiàn)。上述波動(dòng)分析方法是在介質(zhì)為線彈性的條件下進(jìn)行的,為考慮介質(zhì)的非線性性質(zhì),本文利用工程中廣泛應(yīng)用的等效線性化方法將上述線性波動(dòng)有限元分析推廣到非線性領(lǐng)域。為了保證計(jì)算區(qū)內(nèi)波動(dòng)有限元分析的計(jì)算精度,離散單元網(wǎng)格尺寸Δx一般應(yīng)當(dāng)滿足Δx≤(1/6-1/10)λmin(5)Δx≤(1/6?1/10)λmin(5)其中λmin為分析中所關(guān)心頻段內(nèi)諧波的最短波長(zhǎng)。本計(jì)算數(shù)值積分主要采用Newmark-β方法,該方法是無(wú)條件穩(wěn)定的,但在自由場(chǎng)的計(jì)算中,則是依據(jù)差分法的原理來(lái)實(shí)現(xiàn)的,計(jì)算步長(zhǎng)應(yīng)根據(jù)差分法的要求確定。在波動(dòng)有限元分析中,計(jì)算步長(zhǎng)Δt應(yīng)滿足Δt≤α(Δx/c0)min(6)Δt≤α(Δx/c0)min(6)式中c0為場(chǎng)地介質(zhì)的物理波速,一般可取介質(zhì)的剪切波速;α為小于1的系數(shù),根據(jù)實(shí)際計(jì)算情況而定。1.3為相干函數(shù)根據(jù)隨機(jī)過(guò)程理論的基本概念,相干函數(shù)[Coherency]定義為任意2個(gè)隨機(jī)過(guò)程之間的標(biāo)準(zhǔn)化互功率譜密度,即對(duì)于i、j2個(gè)點(diǎn)地震動(dòng),其相應(yīng)相干函數(shù)γij(ω)的嚴(yán)格數(shù)學(xué)定義為γij(ω)={Sij(ω)√Sii(ω)Sjj(ω)0當(dāng)Sii(ω)Sjj(ω)≠0當(dāng)Sii(ω)Sjj(ω)=0(7)γij(ω)=?????Sij(ω)Sii(ω)Sjj(ω)√0當(dāng)Sii(ω)Sjj(ω)≠0當(dāng)Sii(ω)Sjj(ω)=0(7)式中Sij(ω)為2點(diǎn)地震動(dòng)隨機(jī)過(guò)程之間的互功率譜密度;Sii(ω)和Sjj(ω)分別為對(duì)應(yīng)于i、j2個(gè)點(diǎn)的自功率譜密度。相干函數(shù)的另一常用表示方式為γij(ω)=|γij(ω)|exp[iθij(ω)](8)γij(ω)=|γij(ω)|exp[iθij(ω)](8)式中θij(ω)為相干函數(shù)相角;|γij(ω)|稱為遲滯相干函數(shù),用于衡量2個(gè)站點(diǎn)間所記錄運(yùn)動(dòng)過(guò)程線性相關(guān)程度,當(dāng)其中一站點(diǎn)的記錄過(guò)程可由另一站點(diǎn)的記錄經(jīng)某種線性變換得到時(shí),遲滯相干函數(shù)值為1,對(duì)于2個(gè)完全不相關(guān)的過(guò)程,其遲滯相干函數(shù)值為0。當(dāng)?shù)卣饎?dòng)為隨機(jī)過(guò)程時(shí),系統(tǒng)方程為隨機(jī)振動(dòng)方程,本文采用林家浩提出的隨機(jī)振動(dòng)分析的虛擬激勵(lì)方法。在構(gòu)造虛擬激勵(lì)時(shí),根據(jù)波動(dòng)有限元的特點(diǎn),通過(guò)入射波的位移譜來(lái)構(gòu)造虛擬激勵(lì)即虛擬入射波。在入射的地震波為一致激勵(lì)且為直入射情況下,設(shè)入射波位移譜為S(ω),則所構(gòu)造的虛擬入射波為√S(ω)eiωtS(ω)????√eiωt。由確定性結(jié)構(gòu)分析方法得到該虛擬激勵(lì)作用下的場(chǎng)地響應(yīng)矢量y,最后可以寫(xiě)出穩(wěn)態(tài)響應(yīng)譜密度矩陣:Yyy(ω)=y*(ω)yΤ(ω)(9)Yyy(ω)=y?(ω)yT(ω)(9)將計(jì)算得到的波動(dòng)觀測(cè)量的譜密度代入式(7)或(8),即可獲得所求的地震動(dòng)相干函數(shù)。當(dāng)需要考慮入射波的行波效應(yīng)時(shí),入射波位移譜為s(ω)exp(-iωd/v0),其中d,v0分別為任意2點(diǎn)間距離和入射波視波速。設(shè)底部介質(zhì)剪切波速為v,依據(jù)式(10)可確定虛擬剪切波的入射角θ,從而可按文獻(xiàn)的分析方法通過(guò)對(duì)斜入射波的計(jì)算求得問(wèn)題的解答。θ=arcsin(v/v0)(10)θ=arcsin(v/v0)(10)當(dāng)需要考慮土的非線性效應(yīng)時(shí),非線性單元的等效剪應(yīng)變可由下式求得:(γd)eff=0.65pmσγd(11)(γd)eff=0.65pmσγd(11)式中pm為峰值因子,pm的取值與地震動(dòng)的強(qiáng)度、輸入地震動(dòng)譜的類型、土的非線性性質(zhì)均有關(guān)系(文獻(xiàn)),考慮到下面將要計(jì)算問(wèn)題的具體情況,在本文中pm取為2.5;σγd為單元剪應(yīng)變方差,由下式求得:σ2γd=2∫∞0sγ(ω)dω(12)σ2γd=2∫∞0sγ(ω)dω(12)式中sγ(ω)為單元剪應(yīng)變譜密度,可由式(11)的譜密度矩陣及單元的應(yīng)力—應(yīng)變關(guān)系得到。隨機(jī)非線性波動(dòng)分析方法的基本步驟除單元的等效剪應(yīng)變需通過(guò)虛擬激勵(lì)法由上述2式求得外,其他步驟與確定性非線性波動(dòng)分析方法步驟基本一致。依據(jù)波動(dòng)有限元散射問(wèn)題的計(jì)算特點(diǎn),其底部波動(dòng)輸入均采用位移輸入,相應(yīng)的底部波動(dòng)輸入的譜也采用位移譜,當(dāng)?shù)撞坎▌?dòng)輸入的加速度譜已知時(shí),可利用式(13)將其轉(zhuǎn)換為位移譜,然后進(jìn)行波動(dòng)輸入。S¨x¨x(ω)=ω4Sxx(ω)(13)Sx¨x¨(ω)=ω4Sxx(ω)(13)式中S¨x¨x(ω)Sx¨x¨(ω)為入射波加速度自譜密度;Sxx(ω)為入射波位移自譜密度。2分析實(shí)例2.1數(shù)值模擬及結(jié)果分析計(jì)算模型如圖1,土層分2層。表層非線性土層厚30m,計(jì)算物理參數(shù):泊松比為0.3,質(zhì)量密度ρ=2000.0kg/m3。底部硬土層為線性土層,計(jì)算物理參數(shù):泊松比為0.25,質(zhì)量密度ρ=2400.0kg/m3,剪切模量G=6×e8Pa,剪切波速500m/s。表層土的非線性特性以雙曲線模型來(lái)描述:G=Gmax1+|(γd)eff/γR|ξ=|(γd)eff/γR|ξm1+|(γd)eff/γR|(14)G=Gmax1+|(γd)eff/γR|ξ=|(γd)eff/γR|ξm1+|(γd)eff/γR|(14)式中G、ξ、(γd)eff分別為等效剪模、等效阻尼比和等效剪應(yīng)變;ξm和γR為參數(shù),為方便計(jì)算,設(shè)表層土動(dòng)剪切模量的初值G0=Gmax,沿深度一致分布,Gmax=1.2×e8Pa,ξm=0.3,γR=0.0003。底部硬土層內(nèi)入射波的加速度譜采用金井清譜:s(ω)=1+4β2(ω/ω0)2[1-(ω/ω0)2]2+4β2(ω/ω0)2s0(15)s(ω)=1+4β2(ω/ω0)2[1?(ω/ω0)2]2+4β2(ω/ω0)2s0(15)取s0=0.0008,ω0=16.9,β=0.94,相當(dāng)于入射波最大加速度0.06g(硬土層表面最大加速度0.12g)。隨機(jī)波動(dòng)分析時(shí)利用式(13)將上述加速度譜轉(zhuǎn)換成位移譜。人工邊界劃定的計(jì)算范圍200m×44m,單元大小2m×2m,計(jì)算時(shí)步0.0015s。不考慮行波效應(yīng)時(shí),地震波直入射,按隨機(jī)非線性波動(dòng)分析方法求得表面點(diǎn)的加速度譜密度曲線如圖2。為驗(yàn)證本文提出的隨機(jī)波動(dòng)非線性分析方法的正確性,對(duì)上述算例同時(shí)進(jìn)行了數(shù)值模擬。對(duì)輸入地震動(dòng)位移過(guò)程的模擬采用三角級(jí)數(shù)形式:x(t)=Ν∑i=1aicos(ωit+φi)a2i=4.0?S(ω)/ω4?Δω(16)x(t)=∑i=1Naicos(ωit+φi)a2i=4.0?S(ω)/ω4?Δω(16)其中φi為在[0,2π]上均勻分布的隨機(jī)變量;s(ω)為輸入地震動(dòng)加速度過(guò)程譜密度。根據(jù)式(16)生成一定數(shù)量的隨機(jī)地震波樣本,按確定性波動(dòng)非線性分析方法計(jì)算每一地震波樣本輸入下的場(chǎng)地動(dòng)力響應(yīng),并綜合所有結(jié)果計(jì)算得出動(dòng)力響應(yīng)的統(tǒng)計(jì)特征。數(shù)值模擬得到的表面點(diǎn)的位移和加速度的均值及方差過(guò)程如圖3,數(shù)值模擬次數(shù)為500次。由圖3可見(jiàn),表面點(diǎn)的加速度的均值反應(yīng)均趨于0,這與輸入地震動(dòng)為0均值平穩(wěn)隨機(jī)過(guò)程相符;同時(shí)可見(jiàn)隨著時(shí)間的增加,表面點(diǎn)的加速度方差過(guò)程趨于穩(wěn)定。根據(jù)由本文提出的隨機(jī)波動(dòng)非線性分析方法得到的圖2的加速度譜密度s(ω)曲線,按式(12)計(jì)算表面點(diǎn)穩(wěn)態(tài)反應(yīng)的加速度方差,可得σ=0.39m/s2。將此結(jié)果與圖3穩(wěn)態(tài)段的值進(jìn)行比較,可見(jiàn)兩者符合較好,比較結(jié)果說(shuō)明本文提出的隨機(jī)波動(dòng)非線性分析方法是可行的。以下研究不考慮行波效應(yīng)條件下(直入射)場(chǎng)地地震動(dòng)的相干性,場(chǎng)地表面遲滯相干函數(shù)及相角如圖4。由圖4可見(jiàn),表面各點(diǎn)之間的遲滯相干函數(shù)總為1,這與一致激勵(lì)隨機(jī)振動(dòng)基本理論相符;相應(yīng)的相角趨近于0,這與入射波不考慮行波效應(yīng)是一致的。因此對(duì)水平向均勻分布的成層介質(zhì)場(chǎng)地,當(dāng)?shù)卣鸩ù怪比肷鋾r(shí),表面各點(diǎn)之間的遲滯相干函數(shù)總為1,相應(yīng)的相角為0。此外,由圖2可見(jiàn),表面加速度自譜密度具有較窄的頻段,這在一定程度上反映了土介質(zhì)的非線性性質(zhì),計(jì)算同時(shí)表明,在系統(tǒng)的基頻處,表面點(diǎn)加速度譜密度達(dá)到最大值。當(dāng)考慮行波效應(yīng)時(shí),入射波譜具有形式s(ω)·exp(-iωd/v0),假設(shè)入射波視波速v0=1500m/s,底部介質(zhì)的剪切波速v=500m/s,可得入射波與垂直方向夾角為20°。斜入射時(shí)表面點(diǎn)加速度譜密度曲線見(jiàn)圖2,遲滯相干函數(shù)與圖4完全一致,遲滯相干函數(shù)相角變化見(jiàn)圖5。由圖5可見(jiàn),地震波斜入射時(shí),表面各點(diǎn)之間的遲滯相干函數(shù)總為1,相應(yīng)的相角與行波效應(yīng)的相角一致。綜上所述,在一致激勵(lì)情況下,當(dāng)表面介質(zhì)沿水平向均勻分布時(shí),表面各點(diǎn)之間的遲滯相干函數(shù)總為1,相應(yīng)的相角由入射波的行波效應(yīng)決定。應(yīng)該指出,在波動(dòng)有限元分析中,計(jì)算結(jié)果的可靠頻段范圍一般可由式(5)進(jìn)行估計(jì)。設(shè)計(jì)算中關(guān)心的最高頻率為fmax,將式(5)改寫(xiě)為fmax≤(1/6-1/10)(c/Δx)min(5b)fmax≤(1/6?1/10)(c/Δx)min(5b)由上式可見(jiàn),計(jì)算結(jié)果的有意義頻段范圍與單元大小劃分有密切關(guān)系。本例中,表層非線性土層經(jīng)疊代計(jì)算后所得最小剪切模量G=0.4×e8Pa,相應(yīng)剪切波速約為141m/s,可靠頻段的上限約為(1/8)·(141/2)=8.8(Hz)。考慮到場(chǎng)地運(yùn)動(dòng)的高頻成分對(duì)非線性土層的力學(xué)指標(biāo)影響有限,因此按本例單元?jiǎng)澐址绞接?jì)算得到的土層非線性力學(xué)指標(biāo)是合理的,場(chǎng)地地震動(dòng)相干性計(jì)算結(jié)果在有意義的頻段范圍內(nèi)是可靠的。當(dāng)需要在更高的頻段范圍內(nèi)考察場(chǎng)地的地震動(dòng)相干性時(shí),波動(dòng)有限單元網(wǎng)格尚需進(jìn)一步細(xì)劃,這必然會(huì)增加計(jì)算工作量,有關(guān)這一問(wèn)題還需進(jìn)一步研究。2.2計(jì)算方法的穩(wěn)定性分析考慮一側(cè)向介質(zhì)分布不均勻場(chǎng)地(圖6),表層為非線性土層,其下為硬土層,各土層物理力學(xué)指標(biāo)及底部輸入地震動(dòng)譜與算例一相同?？紤]行波效應(yīng),入射波視波速v0=1500m/s,此時(shí)入射波與垂直方向夾角為20°。由于多次透射公式(4)的應(yīng)用條件是在人工邊界附近介質(zhì)沿邊界法向均勻分布,根據(jù)本模型的特點(diǎn),當(dāng)右人工邊界距不規(guī)則地形足夠遠(yuǎn)時(shí),右人工邊界附近非線性介質(zhì)的分布是可以滿足這一條件的。計(jì)算表明為滿足這一條件,由人工邊界劃定的計(jì)算范圍不小于200m×44m,單元大小2m×2m,隨著計(jì)算范圍在水平方向的增大,場(chǎng)地點(diǎn)3、4、5等的反應(yīng)基本不變。這表明只要人工邊界劃分的范圍足夠大,本文的計(jì)算方法是穩(wěn)定的?，F(xiàn)將計(jì)算范圍設(shè)為350m×44m,對(duì)上述計(jì)算范圍進(jìn)行隨機(jī)非線形性波動(dòng)分析,非線性介質(zhì)表面點(diǎn)的加速度譜密度曲線見(jiàn)圖7,圖中同時(shí)繪出非線性介質(zhì)側(cè)向均勻分布場(chǎng)地表面加速度譜密度曲線以資比較。由圖7可見(jiàn)隨著表面點(diǎn)遠(yuǎn)離不規(guī)則地形,其加速度譜密度曲線逐漸接近規(guī)則場(chǎng)地的加速度譜密度曲線。以下研究場(chǎng)地地震動(dòng)的相干性。計(jì)算表明,一致激勵(lì)下表面各點(diǎn)間遲滯相干函數(shù)總為1,與算例一結(jié)果完全相同,故不繪出。表面各點(diǎn)間遲滯相干函數(shù)相角變化情況見(jiàn)圖8、圖9,對(duì)于側(cè)向不均勻分布場(chǎng)地,遲滯相干函數(shù)相角的變化與行波效應(yīng)不一致,與算例一完全不同,這與局部場(chǎng)地條件的影響密切相關(guān)。設(shè)場(chǎng)地表面2個(gè)點(diǎn)i、j的加速度譜密度分別為si(ω)和sj(ω),據(jù)此可得√si(ω)/sj(ω)si(ω