弹性地基梁或板模型被广泛应用于土工加筋路堤的沉降与道路结构受荷载作用下的响应研究中．在这些研究中，地基往往被简化成水平向均匀的Winkler地基[1-3]、多参数地基[4-7]和弹性地基[8-9]．事实上地基不仅沿深度方向存在分层现象，在水平方向上也存在非均匀的情况，如双向增强复合地基[10]和部分机场跑道的非均匀弹性地基[11]． 针对水平向非均匀弹性地基的研究非常有限，大多基于简化的变刚度Winkler地基模型或采用有限元的方法进行分析．文献[12]将复杂地基假设为变刚度Winkler地基，推导出了变刚度Winkler地基梁的弯曲变形解析解；文献[13-15]将复合地基简化为一系列不同刚度的弹簧，对非均匀Winkler地基上梁受静载作用下的变形进行了研究；文献[16]采用微分变换法，对非均匀Winkler地基上变厚度板的自由振动进行了分析．尽管变刚度Winkler地基模型使问题得到简化，但文献[13]的研究表示，距离荷载中心越远，梁的变形结果误差越大，这主要是由于Winkler地基不能考虑地基土的剪切作用．针对这一不足，文献[17]提出半解析有限元法，对水平向非均匀挖填道基受飞机荷载作用下的响应进行了研究．文献[16，18]指出有限元法求解该类问题需要较多的网格和较大的计算量才能达到要求的精度，文献[18]还指出谱方法的优势与应用．谱方法已经被广泛应用于研究弹性地基梁的变形问题[3,5-9]．该方法将解表示为傅里叶级数或利用傅里叶变换的方法进行求解，具有极高的精度，以此为基础的谱元法[18-19]、薄层法[20-21]等方法由于只须在深度方向对梁和地基进行分层计算，因此兼备精度高和计算量小的优势．然而这些研究只考虑了地基沿深度方向上的非均匀性，而忽略了地基沿水平方向上的非均匀性．目前还未有研究利用谱方法计算水平向非均匀弹性地基上梁的弯曲变形问题．本研究利用谱方法探讨水平向非均匀弹性地基上梁的弯曲变形问题．基于欧拉-伯努利梁理论和简化多层弹性地基理论，引入虚拟单元的概念，并考虑地基沿水平方向上的非均匀性，推导出了水平向非均匀弹性地基上梁受静载作用下的控制方程，利用傅里叶级数和伽辽金法相结合的方法，得到了梁和地基变形与应力的数值结果，并与有限元的结果进行了对比，验证了方法的正确性．利用数值结果，进一步分析了地基的水平向非均匀性对梁的弯曲变形的影响．1 控制方程及求解如图1所示，欧拉-伯努利梁和地基完全接触，梁受到的竖向荷载为q(x)(x为X轴坐标)，分布宽度为2l，梁的截面为单位宽度的矩形，厚度和长度分别为hb和2Lb，弹性模量为Eb，抗弯刚度为EbIb，两端为自由边界．地基厚度为hf，底部为基岩，长度为2Lf→∞，沿X方向呈现非均匀性，在图1中以不同的区域标识了出来．10.13245/j.hust.220404.F001图1地基上梁受静载作用的示意图首先对梁进行分析，如图1所示，本研究仅考虑均布荷载的情况，有q(x)=qH [l2-x2]，式中H [∙]为Heaviside函数．假设底部受到地基的反力为p(x)，梁的位移为wb，由于梁的长度为2Lb，其两端为自由边界，因此文献[8-9]引入符合自由边界条件的模态函数进行求解．这种方法的模态函数和计算过程非常复杂．本研究通过引入虚拟梁单元，使得问题大大简化，且该方法有很好的普适性．假设梁分布在整个地基长度，即[-Lf,Lf]上，那么梁的抗弯刚度可以写成：EI=EbIb=Ebhb3/12    (x∈[-Lb,Lb]);0    (x∉[-Lb,Lb]).上式表明，梁在[-Lb,Lb]范围内有刚度，为实体梁，而在该范围外没有刚度，为虚拟梁．此时梁的控制方程可以写成      d2(EId2wb/dx2)dx2=EId4wbdx4+2d(EI)dx∙d3wbdx3+d2(EI)dx2d2wbdx2=q(x)-p(x). (1)由于引入了虚拟梁单元，梁在整个地基长度上抗弯刚度非均匀，因此式(1)中抗弯刚度对x的一阶导数和二阶导数均不为0，即：dEI/dx=EbIbδ(x+Lb)-EbIbδ(x-Lb);d2EIdx2=EbIbd[δ(x+Lb)]dx-EbIbd[δ(x-Lb)]dx.梁的边界条件为wbx=±Lf=∂wb∂xx=±Lf=∂2wb∂x2x=±Lf=∂3wb∂x3x=±Lf=0．梁的位移wb可以写成傅里叶级数的形式，即wb=∑m=-nsnswbmeiαmx;αm=mπLf, (2)式中：2ns+1为能使级数收敛的总项数；wbm为梁位移的傅里叶系数．将式(2)代入式(1)，式(1)两端同乘以e-iαnx，并在-Lf到Lf上积分，可以得到      ∑m=-nsns∫-LfLfEIαm4wbmei(αm-αn)xdx+∑m=-nsns∫-LfLf(-2iαm3)d(EI)dxwbmei(αm-αn)xdx+∑m=-nsns∫-LfLf(-αm2)d2(EI)dx2wbmei(αm-αn)xdx=∫-LfLfq(x)-p(x)e-iαnxdx. (3)δ函数存在以下性质，即      ∫-∞∞djδ(x-x0)dxjf(x)dx=(-1)jdjf(x0)dxj(j=0,1). (4)联立式(3)和式(4)，梁的各振型间的控制方程为      ∑m=-nsns(Abm+Bbm+Cbm)wbm=∫-LfLfq(x)-p(x)e-iαnxdx, (5)式中：Abm=2EbIbαm4sin[(αm-αn)Lb](αm-αn)；Bbm=-4αm3EbIbsin[(αm-αn)Lb]；Cbm=2αm2(αm-αn)EbIbsin[(αm-αn)Lb] .从式(5)可以看出梁的不同级数的系数间存在耦合关系．若Lb=Lf→∞，即为无限长梁的情况，则梁的不同级数系数间相互独立．文献[4-5]提出的简化多层弹性地基模型既考虑了地基土的剪切作用，又通过引入修正参数得到了和有限元一致的结果，可以用于计算弹性地基上梁的弯曲变形．假设地基的弹性模量和泊松比分别为Ef(x)和νf(x)，对应的拉梅常数为λf(x)和Gf(x)，可以表示为：λf(x)=Ef(x)νf(x)[1+νf(x)][1-2νf(x)];Gf(x)=Ef(x)/{2[1+νf(x)]}.假设地基土的竖向位移为wf，地基土的应力可以表示为τxz=FsGf∂wf/∂x;σz=Fn(λf+2Gf)∂wf/∂z,式中：z为Z轴坐标；τxz为剪应力；σz为正应力；Fs和Fn为修正参数[4]， Fs(x)=Fn(x)=0.96-0.49ν(x)+6.47ν(x)2-15.76ν(x)3．根据平衡方程，地基土的控制方程为∂(FsGf∂wf/∂x)∂x+∂(Fn(λf+2Gf)∂wf/∂z)∂z=0．(6)考虑到地基沿水平方向上的非均匀性，式(6)可以写成      ∂(FsGf)∂x∂wf∂x+FsGf∂2wf∂x2+Fn(λf+2Gf)∂2wf/∂z2=0. (7)地基的边界条件可以表示为wfx=±Lf=∂wf∂xx=±Lf=0．竖向位移可以用级数表示为：wf(x,z)=∑m=-nsnswfm(z)eiαmx;αm=mπLf, (8)式中wfm为地基位移的傅里叶系数．同样地，应力也可以写成级数形式为σz=∑m=-nsnsFn(λf+2Gf)dwfmdzeiαmx=∑m=-nsnsσzmeiαmx. (9)同样地，将式(8)代入式(7)，在两边同乘e-iαnx，并在[-Lf,Lf]上积分，可以得到      ∑m=-nsns∫-LfLf∂(FsGf)∂x(iαm)wfmei(αm-αn)xdx+∑m=-nsns∫-LfLfFsGf(-αm2)wfmei(αm-αn)xdx+∑m=-nsns∫-LfLfFn(λf+2Gf)d2wfmdz2ei(αm-αn)xdx=0.从上式可以看出：对于水平向均匀地基，∂(FsGf)/∂x=0，式中的第一项为0；对于水平向非均匀地基，式中的第一项不为0．此时该式可以写成∑m=-nsns(Afm+Bfm)wfm+∑m=-nsnsCfmd2wfmdz2=0，(10)式中：Afm=∫-LfLf∂(FsGf)∂x(iαm)ei(αm-αn)xdx；Bfm=∫-LfLfFsGf(-αm2)ei(αm-αn)xdx；Cfm=∫-LfLfFn(λf+2Gf)ei(αm-αn)xdx．(11)按照文献[20-21]中薄层法的思路，将地基沿深度方向分层，如图2所示，假设共有Nz层，zj表示第j层的深度，wfj，σzj和τxzj分别表示第j层的竖向位移、正应力和剪应力，且假设每一层的位移呈线性变化，即wf=Njwfj+Nj+1wfj+1;wfm=Njwfmj+Nj+1wfmj+1,式中：Nj=zj+1-zzj+1-zj;Nj+1=z-zjzj+1-zj. (12)10.13245/j.hust.220404.F002图2地基上分层示意图将式(12)代入式(10)，根据伽辽金法有：      ∫zjzj+1∑m=-nsns(Afm+Bfm)wfm+∑m=-nsnsCfmd2wfmdz2Njdz=0; (13)      ∫zjzj+1∑m=-nsns(Afm+Bfm)wfm+∑m=-nsnsCfmd2wfmdz2Nj+1dz=0. (14)式(13)和(14)可以分别简化为：      ∑m=-nsns∫zjzj+1[(Afm+Bfm)Nj2-Cfm(dNj/dz)2]∙dz(wfmj)+∫zjzj+1[(Afm+Bfm)NjNj+1-Cfm(dNj/dz)(dNj+1/dz)]dz(wfmj+1)=∑m=-nsnsCfmdwfmjdz;(15)      ∑m=-nsns∫zjzj+1[(Afm+Bfm)NjNj+1-Cfm(dNj/dz)(dNj+1/dz)]dz(wfmj)∫zjzj+1[(Afm+Bfm)(Nj+12)-Cfm(dNj+1/dz)2]dz(wfmj+1)=-∑m=-nsnsCfmdwfmj+1dz. (16)对于上两等式的右边，结合式(9)和(11)，有      ∑m=-nsnsCfmdwfmdz=∫-LfLf∑m=-nsnsFn(λf+2Gf)∙ei(αm-αn)xdxdwfmdz=∫-LfLfσze-iαnxdx.式(15)和式(16)的右边可以分别写成：∑m=-nsnsCfmdwfmjdz=∫-LfLfσzje-iαnxdx；-∑m=-nsnsCfmdwfmj+1dz=-∫-LfLfσzj+1e-iαnxdx．由于欧拉-伯努利梁和地基完全接触，因此梁和地基表层(j=1)之间存在如下关系：p(x)=σz1(x);wb(x)=wf1(x).由于地基底部为基岩，因此底部位移为0，即wfNz+1(x)=0．将Nz每一层的方程按照文献[20]的方法进行组装，并联立式(5)、(15)和(16)，可以得到整个梁-地基系统的控制方程，即Kw=F，(17)式中：K为(2ns+1)Nz×(2ns+1)Nz的矩阵；w和F为(2ns+1)Nz×1的向量．因此，梁的变形可以通过式(2)和式(17)求出，梁的弯矩可以通过级数的形式求出，即Mb=EId2wbdx2=-EI∑m=-nsnsαm2wbmeiαmx．地基土的位移和应力可以通过式(8)、(15)和(17)得出．2 数值验证本研究采用两个算例进行验证．算例1计算水平向均匀地基上的梁的弯曲变形，算例2计算水平向非均匀地基上的梁的弯曲变形．梁的弹性模量为284.4 MPa，厚度为0.4 m[2-3]，长度2Lb=20 m，荷载宽度2l=2Lb=20 m，均布线荷载为q=120 kN/m．利用ABAQUS软件建立二维有限元模型予以验证．算例参数见表1和表2．10.13245/j.hust.220404.T001表1算例1的参数参数弹性模量/MPa泊松比厚度/m地基土1300.2510地基土2300.251010.13245/j.hust.220404.T002表2算例2的参数参数弹性模量/MPa泊松比厚度/m范围/m地基土1300.2510x∈[-∞,0)地基土2600.2510x∈[0,∞]当2Lf=40 m，级数总项数2ns+1=101，沿深度方向分5层时，计算结果收敛．本研究计算的梁的变形结果与ABAQUS的数值模拟结果见图3．10.13245/j.hust.220404.F003图3均布荷载作用下梁的变形结果由图3可见：本研究的计算结果与ABAQUS的计算结果基本符合，仅在梁的两端处有较小的误差，这是由于本研究采用的控制方程仅考虑了梁和地基的竖向位移，而忽略了梁和地基的横向位移．在梁的中心附近，地基的横向位移几乎可以忽略，此时文献[4-5]的简化是非常合理的；而在梁的两端附近，地基土存在较小的横向位移，而对地基的这种假设会产生一定的误差，最大误差约为3%，但总体结果符合得很好，验证了本文方法的正确性．3 参数分析为了进一步分析土体的水平向非均匀性对梁的弯曲变形的影响，本研究对不同水平向非均匀弹性地基上梁的位移、弯矩和竖向应力分布进行了计算．将地基沿水平向分为三部分，如图4所示．其中，地基土1的范围是x∈[-∞,-6] m，地基土2的范围是x∈[-6,6] m，地基土3的范围是x∈[6,∞] m．梁的弹性模量为284.4 MPa，宽度为2Lb=20 m，厚度为hb=0.4 m[2-3]，荷载宽度2l=2Lb=20 m，均布线荷载为q=120 kN/m．所有工况下地基土的泊松比和厚度均相同，即νf1=νf2=νf3=0.25，hf1=hf2=hf3=10 m，不同工况下地基的弹性模量不同，具体参数见表3，工况1为对照组．10.13245/j.hust.220404.F004图4算例示意图10.13245/j.hust.220404.T003表3弹性地基的参数参数/MPa工况1工况2工况3Ef1603060Ef2606030Ef3606060图5～6表示不同工况中梁的弯曲变形结果，图6中Mb为 梁的弯矩．可以看出：工况1的竖向位移最大，工况1和工况3梁的变形和弯矩关于x=0对称，这是由于两个工况中地基的刚度关于x=0对称；不同工况中梁的弯曲变形情况有很大不同，这是由于地基水平向的非均匀性导致的．10.13245/j.hust.220404.F005图5不同工况下梁的变形10.13245/j.hust.220404.F006图6不同工况下梁的弯矩图7表示不同工况中梁底部的应力分布．可以看出：和水平向均匀地基的情况对比，水平向非均匀地基上梁的底部存在应力集中的现象，出现在不同地基土的交界处附近，这主要是由于地基土之间存在刚度差异，而在其他区域的应力约等于均布荷载，这表明应力集中效应的影响范围主要集中在不同地基土的交界处附近，影响范围约为0.5 m，对于本研究的情况，最大应力为均布荷载的1.2倍．10.13245/j.hust.220404.F007图7不同工况下梁底部的应力分布4 结论a. 和有限元法对比，本研究利用谱方法(傅里叶级数)，只在深度方向上进行分层，针对本研究的问题，在深度方向只用分5层，结果就能收敛，兼备了计算量小和结果精度高且易收敛的优点．b. 不同于水平向均匀地基的情况，水平向非均匀弹性地基的不同级数系数之间存在耦合的关系，且地基的水平向非均匀性对地基上梁的弯曲变形有很大影响，在计算中应充分考虑地基沿水平向的非均匀性．c. 水平向非均匀地基上梁的底部存在应力集中现象，出现在不同地基土的交界处附近，这主要是由于不同地基土的刚度差异导致的，对于本文算例，应力集中效应的影响范围约为0.5 m，最大应力为均布荷载的1.2倍．