飽和多孔地基與矩形板動力相互作用的非軸對稱混合邊值問題1)

王立安 趙建昌 楊華中

(蘭州交通大學土木工程學院,蘭州 730070)

引言

地基板(梁)在動力作用下的振動彎曲問題在諸多工程領域有著重要的實踐應用,如建筑物的筏板基礎、公路或市政道路的路面板、鐵路的軌道板以及飛機跑道等都可視為這類問題.以往研究中,通常采用線彈性、彈塑性或黏彈性本構來模擬地基(如Winkler 地基),通過在板?土接觸面上添加彈簧支撐和阻尼黏壺來模擬接觸關系[1-6].眾所周知,土體是一種具有多相組成和骨架孔隙結構的介質,采用彈(塑)性或黏彈性本構來模擬土地基,是將土體視為單相介質,無法真實反映土體的力學性質.尤其在動力作用下,土體中孔隙水的滲透作用以及固流耦合效應對其動力行為起到重要影響[7-10].而且采用添加彈簧支撐和黏壺的方式,使得計算結果完全依賴于地基彈性系數和阻尼比的選取,無法考慮板土接觸的真實情況.采用有限元分析地基板問題時,同樣需要輸入地基彈性系數和阻尼比,而且傳統有限元軟件也沒有提供可模擬滲透和孔隙壓力的本構單元,因此無法實現板土接觸分析的有限元模擬[11-14].何濤[15]基于任意拉格朗日?歐拉(ALE)有限元技術,給出一種數值模擬固流耦合效應的方法,但無法適用于土體.陳少林等[16-17]采用集中質量顯式有限元方法,并結合透射邊界條件給出一種針對海洋地震工程流固耦合問題的統一計算框架.康建宏等[18]基于修正Darcy 模型,分析了多孔介質內黏彈性流體的熱對流問題.Biot[19-20]提出的多孔介質理論為研究土體動力響應提供了理論基礎.Zhang 等[21]給出一種分析飽和土與隧道結構動力相互作用的理論模型.陳少林等[22-23]基于Biot 理論對土?基礎?結構系統在地震作用下的動力響應做了整體分析.

采用接觸分析法求解地基板問題時,關鍵是接觸反力的確定.早期研究中,前蘇聯學者哥爾布諾夫波沙道夫[24]采用雙重冪級數展開給出彈性地基上矩形板的完全解析解,但這種解法形式太過復雜,并沒有在實踐中得到應用.王春玲等[25-26]采用雙重三角級數展開的方法,分別給出了彈性地基和飽和多孔地基上矩形板動力彎曲的解析解.這種方法與波沙道夫的方法相似,形式復雜存在收斂問題,而且無法求解復雜邊界條件的情形.一些學者將矩形板剛性化,假定接觸反力均勻分布,從而使問題簡化為平面應變問題進行求解,對地基上剛性基礎的振動問題做了解析研究[27-29].還有一些學者將矩形板視為無限大使問題簡化為單一邊界問題從而進行解析求解.孫璐等[30]給出彈性半空間上無限板的解析解;蔡袁強等[31-32]給出飽和多孔地基上無限大板的動力彎曲解;張春麗等[33]給出各向異性彈性地基上無限大板的解析解.金波等[34-35]、陳龍珠等[36]利用Abel 變換給出了軸對稱混合邊值問題的求解方法,并由此得出飽和多孔地基上圓形板振動彎曲的解析解.應注意,對于圓形地基板的軸對稱問題,放在圓柱坐標系下事實上是二維問題,其一般性將受到限制.此外,郭樹起[37]也給出了適用于軸對稱問題的邊界積分法.而對應于矩形地基板的三維非軸對稱混合邊值問題的完全解析解,目前仍未見相關文獻.

基于以上分析,本文考慮地基土的固流耦合效應和孔隙水的滲透作用,基于Kirchhoff 理論和Biot 多孔介質理論建立矩形板與飽和多孔地基動力相互作用的計算模型.并且考慮了一種混合透水邊界,即板土接觸面為不透水界面,其余表面為透水界面.針對接觸面和非接觸面的混合邊值問題,通過雙重Fourier變換構造出兩對二維對偶積分方程,并利用Hilbert正交化定理和Schmidt 法對二維對偶積分方程進行求解,最終推導出板土系統在動力作用下的位移和應力解析式.通過數值計算,對飽和多孔地基上矩形板的振動響應和參數影響做了分析和討論.將計算結果與彈性地基上矩形板的振動響應進行對比分析,總結出相關結論.

1 計算模型及基本方程求解

如圖1 所示,長、寬分別為2la和2lb的矩形彈性板支撐在飽和多孔半空間表面.板長方向為x軸,板寬方向為y軸,z軸指向地基.板中心作用有長為2a寬為2b的簡諧分布載荷.圖1(b)表示板?土接觸面為不透水界面,其余地表為透水界面.當板?地基系統受到簡諧載荷激勵時,系統中各場量存在如下關系

圖1 模型示意圖Fig.1 Diagram of calculation model

式中,wb為矩形板撓度,F為接觸反力,u為地基土骨架位移,w為孔隙水位移,σ 和τ 分別為土骨架正應力和切應力,p為孔隙水壓力.為對應物理量的幅值.eiωt為諧和因子,ω 為角頻率,t為時間,i 為虛數單位.

1.1 矩形板基本方程求解

根據Kirchhoff 理論,矩形板撓曲方程表示為

式中,已根據式(1)將諧和因子eiωt約去;D為矩形板的撓曲剛度,?4=?4/(?x4)+2?4/(?x2?y2)+?4/(?y4);H()為Heaviside 階躍函數.對x,y坐標引入以下雙重Fourier 變換

將式(3)用于式(2),解得

1.2 飽和多孔地基控制方程及求解

考慮孔隙流體壓縮性和流體與骨架的相對運動,并忽略體力作用,Biot 飽和多孔介質波動方程寫為[19-23]

本構方程

滲流連續方程

式中,θ=?ux/(?x)+?uy/(?y)+?uz/(?z)為固體骨架的體應變;λ 和G為Lamb 常數;ρ=nρf+(1 ?n)ρs為土體總質量密度,其中ρs和ρf分別為固體顆粒和孔隙水的質量密度;n為孔隙率;kd為滲透系數;α 為表征固體顆粒壓縮性的常數,當考慮顆粒不可壓縮時α=1.

對式(5a)～式(5c)施以?(5a)/(?x)+?(5b)/(?y)+?(5c)/?z運算,并考慮式(7),得到

對式(5d)～式(5f)施以?(5d)/(?x)+?(5e)/(?y)+?(5f)/(?z)運算,得到

將式(8)代入式(9),得到

對式(10)進行雙重Fourier 變換后得到

求解常微分方程(11),得到

式中,r1=,取Re(r1)>0,以滿足波的輻射條件,Re表示取實部;φ=ρfω2/n+iωρfg/kd;M0=(ρω2?φ)/(λ+2G).

對式(9)做Fourier 變換,再將式(12)代入,即可求出孔隙壓力的變換域解

式中,r2=,M1=(ρfω2?φ)/M0.

由式(4d)～式(4f)可得到

將式(14)分別代入式(5a)～式(5c),消去wx,wy和wz后進行雙重Fourier 變換,得到

式中,M2=λ +G+M1M3,M3=ρfω2/? ?1,M4=,?=ρfω2/n?iωρfg/kd.

求解二階常微分方程組(15),得出位移解

式中,r3=,A1～A4為積分常量,由邊界條件確定;對本構方程(6)做雙重Fourier 變換,并將位移解式(16a)～式(16c)代入,得到變換域的應力解

1.3 混合邊值問題求解

考慮矩形板與半空間表面光滑接觸,且不會發生脫離;并考慮板土接觸面為不透水邊界,其余表面為透水邊界.將接觸面區域(|x|≤la,|y|≤lb,z=0)記為?1,接觸面以外區域(|x|>la,|y|>lb,z=0)記為?2,則半空間表面z=0 處的邊界條件寫為

利用式(13)、式(17a)、式(18a)得出

其中,N0及N11～N33見附錄.將式(19)代入式(4)、式(14)、式(16)中并進行Fourier 變換,得到

式中,h,Δ,E1,E2,H1,H2均為ξ 和η 的已函數,形式見附錄.對式(20)做Fourier 逆變換,而后代入式(18a)和式(18b),整理出兩對二維對偶積分方程

為求解對偶積分方程(21),根據Hilbert 正交化定理[38],利用Jacobi 正交多項式將積分核和展為如下級數形式[39]

式中,cmn和dmn為未知系數,為Jacobi 多項式.對式(22)進行雙重Fourier 變換得到[40]

式中

其中,Γ(·),Jn(·)分別為Gamma 函數和Bessel 函數.將式(23b)代入式(21),發現自動滿足,將式(23b)代入式(21)后,則得到

然后利用式(4)、式(16)、式(17)可分別得到矩形板和地基的位移和應力解.

2 收斂性和正確性驗證

首先進行收斂性分析,以確定式(24)的級數項數m和n.為了便于比較,矩形板和地基的計算參數按文獻[20] 取值.取矩形板邊長2la=2lb=4 m,板厚為0.2 m;板的泊松比為0.167,彈性模量Eb=34 300 MPa;地基泊松比為0.4,彈性模量Es=343 MPa;板中心作用均布載荷q=0.98 MPa.本文考慮的地基為兩相飽和多孔地基,所以還需補充土體孔隙率n=0.6,固體骨架密度ρs=2700 kg/m3、孔隙水密度ρf=1000 kg/m3,滲透系數kd=10?5m/s.通過MATLAB 編程,對本文推導結果進行數值計算.圖2給出了矩形板中心點撓度幅值與網格數目m×n的關系曲線,圖中反映出本方法收斂速度很快,只需取m×n=16×16 個網格點足以得到收斂結果.

圖2 結果收斂性試算Fig.2 Verification of convergence

若取孔隙水密度ρf=0,代表地基中無水,則本文的地基模型退化為單一彈性地基,可與文獻[20]的結果進行對應比較.表1 列出了本文退化解與文獻解的對比,從表1 的對比結果可知,本文退化解與文獻[20] 很好吻合,由此驗證了本文推導方法和數值計算程序的正確性.

表1 板中心撓度幅值及彎矩與文獻[20]對比Table 1 Comparison of calculation results

3 數值結果和討論

利用前文確定的計算參數和程序,對本文模型進行數值計算圖3 給出了部分計算結果.圖4 為不同載荷激勵頻率下板中心位置的撓度時程曲線,為了進行對比分析,在圖4(a)中繪出了對應于彈性地基的計算結果.對比圖4(a)和圖4(b)可知,對于彈性地基,載荷激勵頻率的變化會造成板振動頻率的變化,但不會改變其振動幅值.而對于兩相飽和多孔地基,激勵頻率的變化不僅引起板振動頻率的變化,其振動幅值也發生大的變化.出現這一現象的機理可理解為,在載荷作用下,飽和多孔地基中固體骨架發生壓縮變形,孔隙受到擠壓,孔隙壓力進而升高,而孔隙壓力會阻礙固體骨架進一步變形.若孔隙中的水能及時滲透,則孔隙壓力將得以消散,這種阻礙效應將不存在.試想一下,當振動頻率較高時,孔隙水來不及滲透,則孔隙壓力持續升高,這種阻礙效應將增大,地基將難以發生變形,地基剛性從而增大,板的撓曲幅值將隨之減小.

圖3 數值計算結果(ω=10 rad/s,t=0 s)Fig.3 Results of numerical computation with(ω=10 rad/s,t=0 s)

圖4 矩形板中心點撓度時程曲線Fig.4 Variation of deflection with time at the center of rectangular plate

圖5 撓度幅值與激勵頻率關系圖Fig.5 Deflection amplitude relative to frequency

為進一步考察飽和多孔地基上矩形板的振動撓曲與激勵頻率的關系,圖5 給出激勵頻率與撓度幅值的關系曲線.分析圖5 可知,隨著激勵頻率增大,板的撓度幅值隨之減小.當載荷激勵頻率達到極限時,板的撓度幅值收斂于一固定值.圖中還反映出,撓度幅值的收斂速度與地基滲透系數有密切關系.為考察撓度幅值的最終收斂值,嘗試取地基滲透系數kd→0.計算發現,激勵頻率達到極限時,板的撓度幅值與地基滲透系數趨于0 時的結果近似相等(即wb|ω→∞≈wb|kd→0).這說明,當載荷激勵頻率達到極限時,地基土體相當于一個固流耦合的封閉系統(孔隙水不發生滲透).

圖6(a)和圖6(b)分別為矩形板中心點下0 m 和2 m 處地基土中孔壓隨時間的變化曲線.圖7(a)和圖7(b)為矩形板中心點下土體中固體骨架位移幅值U和孔隙壓力幅值P沿深度的分布曲線.圖中可以看出,隨著載荷激勵頻率增大,地基中固體骨架的水平方向和豎直方向位移幅值都隨之減小,位移響應深度也隨即減小.但孔隙壓力幅值卻隨激勵頻率的增大而增大,這一結果也印證了前文的機理解釋.

圖6 板中心點下地基中空隙壓力時程曲線Fig.6 Variation of pores pressure with time below the center point of rectangular plate

圖6 板中心點下地基中空隙壓力時程曲線(續)Fig.6 Variation of pores pressure with time below the center point of rectangular plate(continued)

圖7 地基中變形和孔隙壓力Fig.7 Displacement and pore pressure of foundation

圖8 為激勵頻率ω=30 rad/s 時,土?板模量比(Es/Eb)與板的撓度幅值的關系.圖中可看出,Es/Eb取值越小,地基土相對于板越松軟,板的變形越大.圖7(b)比較了相同模量比取值下彈性地基與飽和多孔地基上的板撓曲變形,進一步說明了飽和多孔地基中孔隙壓力對地基剛性的貢獻.

圖8 模量比(Es/Eb)對地基板撓度的影響Fig.8 Influence of modulus ratio on plate deflection

4 結論

本文建立了具有混合透水邊界的板土接觸模型,利用Hilbert 定理和Schmidt 法解決了非軸對稱混合邊值問題的求解,從而得出飽和多孔地基上矩形板振動彎曲的解析解.通過數值計算,發現飽和多孔地基中由于孔隙壓力的作用,地基剛性將隨激勵頻率的增大而提高.而地基板的振動變形將隨載荷激勵頻率的增大而減小.本文提出的非軸對稱混合邊值問題的解析法,具有一般性,可廣泛應用于復雜接觸問題和多場耦合問題的求解.

附 錄