陳少林 孫 杰 柯小飛
(南京航空航天大學土木與機場工程系,南京 210016)
隨著我國海洋戰(zhàn)略的實施,將要建造大量近海工程結構,如何保證近海工程結構在服役期的地震安全性是目前面臨的重要科技問題.要保證近海工程結構的地震安全性,首先要對其在潛在地震作用下的反應進行準確的分析預測[1-9].目前,受水下振動臺設備的限制,采用試驗方法對海工結構的地震反應分析還很少.而在數(shù)值分析方面,其屬于半無限域地震波散射問題(如圖1 所示),涉及到海域場地自由場分析(為散射問題提供輸入)、海水-飽和海床-結構間的相互耦合分析,以及半無限海域場地的人工邊界條件,問題十分復雜.
圖1 海洋工程結構地震反應分析示意圖Fig.1 Seismic response analysis of marine engineering structure
其中,海水、飽和海床、結構間的流固耦合問題,若采用常規(guī)的流固耦合求解思路,須分別采用流體聲波方程、飽和多孔介質方程和固體波動方程的求解器,各求解器在每一時步進行界面耦合,十分不便,目前還沒有成熟方法和軟件能高效地解決此問題.因此,發(fā)展一種海洋地震工程中海水-飽和海床-結構耦合分析的高效方法,考慮海床的非線性,揭示強震作用下海床液化失穩(wěn)機理和海工結構的災害模式十分必要.
對于海域場地非線性地震響應分析,可以采用等效線性化分析方法[10-11]或時域非線性分析方法[12-18].一般認為,等效線性化方法適合于弱非線性,對于強非線性,宜采用時域非線性分析方法.當采用時域非線性方法分析場地地震響應時,需要先進行水平成層場地的自由場分析,作為邊界的輸入.目前,一般采用傳遞矩陣方法[19-21]進行自由場分析,適合于線性或等效線性化情形,可以考慮地震波斜入射;也可采用有限元數(shù)值分析方法[22-24],根據(jù)Snell定律,將二維自由場分析一維化,采用數(shù)值方法實現(xiàn)一維分析,再根據(jù)視速度將一維波場外推至二維、三維,這種外推只適合于波速恒定的線性情形.因此,若采用數(shù)值方法求解非線性自由場,只能獲得一維自由場,即垂直入射情形.在考慮飽和海床的非線性時,理論上,自由場計算也需采用與內域計算同樣的本構,若與內域分析的本構不一致,對結果的影響值得討論.
本文在海域場地自由場分析[20]和海洋地震工程流固耦合問題統(tǒng)一計算框架[25-26]的基礎上,首先考慮飽和海床的非線性,分析海域場地的非線性反應特征,并討論了采用線性自由場和非線性自由場對結果的影響.進一步分析了海域場地和海水-海床-結構體系在地震作用下的動力響應,對比了線性海床和非線性海床情形時響應的差異,為海洋工程結構的抗震設計提供理論基礎和分析手段.
如圖2 所示,從場地中截取有限土體,采用八結點三維實體單元對其進行空間離散,整個計算區(qū)域的有限元節(jié)點可以分為兩類:內節(jié)點與人工邊界節(jié)點.內節(jié)點又分為一般內節(jié)點(包括自由表面節(jié)點)與不同介質交界面上的界面節(jié)點.場地介質為飽和多孔介質,海水、基巖均可視為飽和多孔介質的特殊情形.
圖2 場地-結構相互作用分析模型示意圖Fig.2 Site-structure interaction analysis model
對于飽和多孔介質,其固相平衡方程[27-28]
液相平衡方程
相容方程(考慮初始孔壓和初始體應變?yōu)榱銜r)
其中,Ls,Lw為微分算子矩陣,σ′為有效應力矢量,τ為平均孔壓,P為孔隙水壓,為固、液相的速度,為固、液相的加速度,ρs,ρw分別為固、液相的密度,β 為孔隙率,b=β2μ0/k0,k0為流體滲透系數(shù),μ0為動黏度系數(shù),Ew為流體的體變模量,es,ew分別為固相和液相的體應變.
以上方程當孔隙比β=1 時,可退化為理想流體方程;當孔隙比β=0 時,可退化為彈性波動方程.因此,流體、彈性固體、飽和多孔介質可統(tǒng)一由廣義飽和多孔介質(0 ≤β ≤1)方程描述,僅是材料參數(shù)不同.因此,海域場地中地震波傳播問題可由廣義飽和多孔介質方程統(tǒng)一描述.
對方程(1)和(2)采用伽遼金方法離散,并考慮相容方程(3)及邊界條件,可得任一結點i的運動平衡方程為[26]
其中,Msi,Mwi分別為固、液相質量;分別為固、液相本構力,本文采用Davidenkov 模型考慮固相的非線性,具體見1.4 節(jié);分別為固、液相黏性阻力;分別為固、液相邊界力,當節(jié)點i為一般內節(jié)點時,界面力為零.
對式(4a),(4b)進行時步積分,可得到一般內節(jié)點i的固、液相位移遞推公式
對于界面節(jié)點,如圖2 中的j(k),根據(jù)文獻[26],采用隔離體概念,對式(4a),(4b)進行時步積分得
其中,i=j,k
由式(6)、式(7),以及界面連續(xù)條件,經(jīng)推導,可得
其中
為了有效地模擬外行波穿過人工邊界的運動狀態(tài),我們采用多次透射公式[29]
該局部人工邊界條件具有普適性,與具體波動方程無關,可直接用于飽和多孔介質中的波動[30].分別將式(11)應用于廣義飽和多孔介質的固相和液相位移,可得到邊界點在p+1 時刻的散射波場位移,再疊加上邊界點的自由場位移,即可得到邊界點的固、液相總位移.
在強震作用下,土體將進入非線性,本文采用三參數(shù)的Davidenkov 本構模型,通過對飽和土體固相剪切模量實時修正,可以較為有效地描述飽和土體的非線性動力特性[14-15].
由Martin 和Seed[18]提出的Davidenkov 本構模型骨架曲線的表達式為
式中,τ,γ 分別為剪應力、剪應變;Gmax為初始剪切模量;A,B,γ0可由土體動三軸實驗數(shù)據(jù)進行擬合得到.
經(jīng)Masing 法則修正后滯回曲線如圖3 所示,依據(jù)Pyke 提出的“n倍法”,應力-應變滯回曲線服從下述關系
分別對式(13)、式(15)求導,可得到初始骨架曲線段、滯回曲線段的時變切線剪切模量
圖3 不規(guī)則加卸載準則修正的Davidenkov 本構模型曲線Fig.3 Davidenkov mode modified by irregular loading and unloading criterion
式中參數(shù)(2nγ0)2B根據(jù)當前拐點及歷史上最大(小)點確定.上述本構模型是依據(jù)實驗室得到的主空間上的主剪切應力-主剪切應變曲線而建立的,通常假定該本構模型同樣適合于復雜三維情形的主剪切應力-主剪切應變關系,因此將其用于三維計算分析時,需根據(jù)應變分量,求出主剪切應變,再根據(jù)式(16)和式(17)更新剪切模量.
考慮如圖4 所示的三維海水-飽和土-基巖海域場地模型,整體尺寸為62 m×62 m×60 m,單元尺寸為1 m×1 m×1 m,各層介質參數(shù)如表1 所示.其中,飽和土層采用三參數(shù)的Davidenkov 非線性本構模型,參數(shù)A=1.1,B=0.2,γ0=0.000 4;海水層和基巖層考慮為線性.考慮SV 波垂直入射,輸入如圖5所示的單位脈沖,步距為0.000 2 s,步數(shù)為8000 步.
圖4 海水-飽和土-基巖場地模型示意圖Fig.4 Water-saturated soil-bedrock site model
表1 模型參數(shù)Table 1 Parameters of the model
圖5 脈沖時程圖Fig.5 Pulse time history
分別采用如下兩種方案計算自由場:(1)方案I采用傳遞矩陣方法,得到線性自由場;(2)方案II 采用一維有限元方法,飽和土層采用Davidenkov 非線性本構模型,參數(shù)與三維模型相同,底邊界采用透射邊界,得到非線性自由場.將自由場作為輸入,并在5個邊界面(左、后、右、前、底)上施加透射邊界,進行三維非線性海域場地分析.
理論上,本模型為水平成層場地,SV 波從底部垂直入射,忽略海水黏性,因此SV 波不能在海水層傳播,海水層X向位移應為0.從圖6(a)、圖6(b)中可以看到,方案II 海水層表面A點X向位移為0,與理論相符;而方案I 中,A點X向位移明顯,幅值量級與輸入脈沖波相同,與理論不符.
圖7 為方案II 中,一維與三維模型中飽和土層和基巖層各監(jiān)控點(見圖4)的響應,從圖中可以看出,一維模型和三維模型的響應完全一致,驗證了非線性自由場輸入的可靠性.
圖6 方案I、II 海水層監(jiān)控點X 向位移對比Fig.6 Comparison of X-direction displacement of water layer monitoring points in scheme I and II
圖6 方案I、II 海水層監(jiān)控點X 向位移對比(續(xù))Fig.6 Comparison of X-direction displacement of water layer monitoring points in scheme I and II(continued)
圖7 方案II 一維模型與三維模型X 向位移對比Fig.7 Comparison of X-direction displacement between 1D model and 3D model in scheme II
采用與2.1 節(jié)中相同的模型、脈沖輸入,分別考慮飽和土層線性和非線性情形,分析飽和土層非線性對海域場地響應的影響.
對于線性情形,SV 波在場地中的傳播過程如圖8 所示,根據(jù)基巖層剪切波速c3和土層剪切波速c2估算各波峰到時,Δt3=H3/c3=0.075 0 s,Δt2=H2/c2=0.084 4 s.
圖8 SV 波傳播過程示意圖Fig.8 Schematic diagram of SV wave propagation
根據(jù)文獻[31-32],得SV 波垂直入射時,略去高階項后,基巖層-飽和土層分界面上行波反射系數(shù)
其中波阻抗比ψ=ρ2c2/(ρ3c3),ρ2=ρs2(1 -β2),ρ3=ρs3,得R1=0.547.
上行波透射系數(shù)
得基巖分界面處透射系數(shù)為T1=1.547;基巖層-飽和土層分界面下行波反射系數(shù)
其中波阻抗比ψ′=ρ3c3/(ρ2c2),得R2=-0.547;下行波透射系數(shù)
得基巖分界面處透射系數(shù)為T2=0.453;類似地,得飽和土層-海水層分界面上行波反射系數(shù)為1.
表2 與表3 對波峰到時及幅值的實際值與理論值進行對比,結果誤差較小,證明了線性情形模擬準確.
考慮飽和海床為非線性時,當SV 波傳至飽和土層時,土骨架進入非線性狀態(tài),剪切模量G2減小,則剪切波速c2減小,因此波阻抗比ψ 減小,且0 <ψ <1,根據(jù)式(18)、式(19),反射系數(shù)R1增大,透射系數(shù)T1也增大,因此圖9 中波峰2′,3′,4′峰值比線性情形的要大;圖9(a)中波峰5′峰值本應隨波峰3′峰值增大而增大,但由于飽和土層非線性阻尼效應的影響,峰值反而減小;同時,由于剪切波速c2減小,圖9 中波峰5′,6′,7′,8′到時相比線性情形有延遲;另一方面,剪切波速c2減小,則ψ′增大,且ψ′>1,根據(jù)式(21),基巖層-飽和土層交界面下行波透射系數(shù)T2減小,再加上飽和土層非線性的阻尼效應的影響,波峰6′,7′,8′峰值減小明顯.
此外,對于飽和土層液相,黏性力與慣性力相比,占主導作用.因此,圖9(b)和圖9(d)中液相位移時程表現(xiàn)為明顯的擴散過程.
考慮如圖10 所示的三維場地-結構模型,場地尺寸為62 m×62 m×60 m,樁總長60 m,海水以上部分長10 m,截面尺寸為4 m×4 m,上部承臺尺寸為6 m×6 m×6 m.單元尺寸為1 m×1 m×1 m.飽和土層分別按線性、非線性兩種情形計算,場地、樁和承臺參數(shù)見表1.輸入圖5 所示的的SV 脈沖波.
表2 波峰到時Table 2 Arrival time of wave crest
表3 波峰幅值Table 3 Amplitude of wave crest
圖9 線性和非線性情形各監(jiān)控點X 向位移對比Fig.9 Comparison of X-direction displacement of monitoring points in linear and nonlinear cases
圖9 線性和非線性情形各監(jiān)控點X 向位移對比(續(xù))Fig.9 Comparison of X-direction displacement of monitoring points in linear and nonlinear cases(continued)
圖10 三維場地-結構模型示意圖Fig.10 3D site-structure model
由于樁與結構的存在,樁的運動會弓起海水的運動,從圖11(a)和圖11(b)中可以看到,海水層存在一定位移.
從圖11(c)~圖11(f)可以看出,相比線性情形,非線性情形飽和土層位移幅值要小,衰減更快;非線性情形基巖層位移主要是初始的入射波及其在基巖層-飽和土層界面的反射波,來自飽和土層的透射波幅值很小,原因與2.2 節(jié)中的分析相同.
從圖11(g)~圖11(j)可以看出,由于飽和土層位移影響,非線性情形承臺頂部、樁上部位移幅值要比線性情形小.從圖11(l)可以看出,由于基巖層位移影響,非線性情形樁底部0.3 s 后位移較小.
根據(jù)樁單元節(jié)點的位移,可以得到不同高度處樁水平截面上的剪力和彎矩,如圖12 所示.從圖中可看出,非線性飽和海床與線性飽和海床情形相比,沿樁身的最大剪力和最大彎矩整體上要小.無論是線性還是非線性情形,樁頭的彎矩(-10 m 處)均較大,且在土層界面(10 m 和30 m 處)附近彎矩為極大值,樁端附近彎矩也較大.
本算例中飽和土層厚20 m,在飽和土層中間處水平截取一個平面,監(jiān)控該截面不同時刻的孔壓.圖13 為飽和土層為非線性情形下,SV 波垂直入射時,不同時刻該截面的孔壓云圖.
考慮SV 波沿X–Z平面垂直入射時,主要反應為X方向,由于模型對稱,則位移關于樁基X軸、Y軸對稱分布.孔壓值主要受速度散度變化的影響,孔壓的分布關于樁基X軸對稱,關于Y軸反對稱分布,如圖13 所示.經(jīng)計算地震波從基巖層入射,響應峰值0.25 s 左右到達監(jiān)控面處.下圖中可以看出0.05~0.25 s 階段,由于樁-土位移差異,樁周圍孔壓逐漸集中;0.25 s 左右樁周圍孔壓到達峰值;0.25~1.6 s高壓區(qū)逐漸由樁周圍向四周擴散,孔壓逐漸減小,整個截面孔壓趨于均勻.
圖11 監(jiān)測點X 向位移對比Fig.11 Comparison of X-direction displacement of monitoring points in linear and nonlinear cases
圖11 監(jiān)測點X 向位移對比(續(xù))Fig.11 Comparison of X-direction displacement of monitoring points in linear and nonlinear cases(continued)
圖12 樁截面內力最大值隨深度變化曲線Fig.12 Maximum of internal force with depth
圖13 不同時刻監(jiān)控截面孔壓云圖(單位為Pa)Fig.13 Distribution of pore water pressure on monitoring section(unit:Pa)
圖13 不同時刻監(jiān)控截面孔壓云圖(單位為Pa)(續(xù))Fig.13 Distribution of pore water pressure on monitoring section(unit:Pa)(continued)
本文發(fā)展了一套海域場地和海水-飽和海床-結構體系地震響應的分析方法,內域節(jié)點采用海水、飽和海床、基巖流固耦合分析的統(tǒng)一計算框架,采用透射邊界模擬無限域的影響,分析了SV 波垂直入射時,海域場地和海水-飽和海床-結構體系的響應.討論了海域場地非線性響應分析時,自由場對結果的影響,分析了線性飽和海床和非線性飽和海床情形時響應的差異.得到如下結論:
(1)就本文算例而言,計算海域場地非線性地震響應時,線性自由場將弓起較大的誤差,應采用非線性自由場,且自由場分析的本構應與海域場地分析的本構一致.
(2)由于飽和土層進入非線性,模量減小,由基巖入射到飽和土層的反射系數(shù)和透射系數(shù)均增大,基巖的響應比線性情形時要大.透射系數(shù)增大,飽和土層的反應理應增大,但由于非線性時阻尼增大,導致飽和土層的非線性響應比線性情形時要小,且持時要短.
(3)就本文算例而言,非線性飽和海床情形沿樁身的最大彎矩和剪力要比線性飽和海床情形的小,非線性飽和海床弓起的結構破壞可能主要體現(xiàn)在海床液化弓起的基礎失效方面.