饱和多孔地基与矩形板动力相互作用的非轴对称混合边值问题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 法解决了非轴对称混合边值问题的求解,从而得出饱和多孔地基上矩形板振动弯曲的解析解.通过数值计算,发现饱和多孔地基中由于孔隙压力的作用,地基刚性将随激励频率的增大而提高.而地基板的振动变形将随载荷激励频率的增大而减小.本文提出的非轴对称混合边值问题的解析法,具有一般性,可广泛应用于复杂接触问题和多场耦合问题的求解.

附 录