基于无单元galerkin法的超长嵌岩灌注桩屈曲稳定性分析

基于无单元galerkin法的超长嵌岩灌注桩屈曲稳定性分析

洞庭湖区的顶层很脆弱很深。为了满足上部桥梁结构对基本承受和变形的要求，主要采用桩基，通常是高承台式的长石桩。而这类桩基的承载力取决于桩顶容许变形及桩身压屈稳定性。事实上,当桩周上覆土体软弱或易振动液化、桩顶自由长度较大(高承台)及桩身较长且嵌岩时,桩身可能如细长杆产生突发且后果严重的屈曲破坏。但直至20世纪20年代,有关桩的屈曲稳定问题才引起国内外学者重视,如Lee(1968)的室内模型钢铝桩试验结果证明了软土中桩身屈曲破坏的可能性,Reddy(1970)等提出的屈曲分析能量法解答。国内70年代胡人礼等首先提出了一组桩身计算长度经验公式(现仍被采用)。至90年代,随着超长桩的逐渐使用,桩身压屈问题受到进一步重视,不少学者开展了更深入的研究,如赵明华提出的桩身稳定长度计算公式及所完成的室内模型试验等。然而,已有方法几乎无一例外是基于第一类稳定问题的特征值法:视桩身为理想轴压杆,不考虑桩-土共同作用,将桩身失稳临界荷载归结为稳定特征方程的求解,所获解答为第二类问题(由于桩身初始缺陷、桩土接触与材料非线性等,桩身屈曲稳定应属第二类)的上限,用于工程偏于不安全,故有必要按第二类问题采用数值法分析桩身屈曲稳定性。作为数值分析的有力工具,有限元法已在工程领域取得巨大成就。但对于基桩屈曲分析中存在的刚度悬殊、滑移与开裂等非变形协调问题,网格的存在妨碍处理与原始网格线不一致的不连续性和大变形,需不断进行网格重构,导致计算结果精度降低甚至不收敛,有必要寻求脱离于单元的数值法。20世纪70年代以来迅速发展的无单元法(meshlessmethod)就是在此需求下产生的,其共同点是:基于点近似,无需划分单元,具有良好的动态追踪裂纹或不连续面扩展的能力,其中尤以无单元伽辽金法(ElementFreeGalerkinMethod,EFGM)具有求解稳定、精度高、并可消除体积闭锁等优点而应用广泛,故本文引入该法进行桩身屈曲分析。1点x支持域+点x-1无单元法基于分析域上某一局部范围内的离散场点采用移动最小二乘(MovingLeastSquare,MLS)法插值后得到近似场函数。设u(x)为分析域内的位移场函数,域内某点x(对二维问题,x=(x,y)T,u=(u,v)T)处的近似值uh(x)可由MLS法表示为:其中:a(x)为与坐标x有关的系数向量;PT(x)为m阶完备多项式基函数向量,对二维问题表示为:系数向量a(x)(=(a1(x),a2(x),...,am(x))T)可由点x支持域(supportdomain)内n个节点(x1,x2,...,xn)的场值求得。所谓点x的支持域,即局部范围内用来求得该点近似函数值的若干节点数,常用的形式有圆形与矩形等,如图1所示。设点x的支持域内有n个节点,其位移分别为u1,u2,...,un,则可由MLS法即式(1)得该节点的局部近似函数值为:为确定a(x),先构建带权重的L2范数J(x):式中,w(x-xI)为非负的具有紧支性的连续可导型权函数,即支持域内各节点的权函数在点x处的权值,常用形式有指数型、样条型等。然后,系数向量a(x)可由式(4)的加权残量最小原则求解,即由∂J/∂a=0可得:将式(5)代入式(1),并整理后可得:式中,φI(x)、ΦT(x)分别为形函数和形函数矩阵。采用随空间坐标变化的移动权函数w(x)及系数向量a(x),是MLS区别于其他最小二乘法的关键所在。若所有权函数所给权值为常量,则式(1)即为标准的非移动最小二乘近似,且a(x)也成为常系数向量。但MLS近似不一定精确通过插值节点,构造的形函数不具备Kroneckerδ函数性质,导致边界条件(essentialboundaryconditions)的施加不如有限元法直接方便,需采用其他辅助法实现。2桩底板的非变形协调材料的不连续处理2.1适用条件分析目前,国内外学者提出了不少接触面本构模型[11―12],各有其特点与适用性,且更多的是剪应力τs与剪切位移s间的关系。经综合考虑,本文采用非线性弹性-理想塑性的桩土接触面本构关系,如图2所示。2.2节点位移向量考虑一般情况,建立桩土位移不连续界面元局部坐标系及其方向如图3所示,其中s轴与整体坐标系x轴夹角为θ,以u′,v′和u,v分别表示局部与整体坐标系下节点的位移,并规定法向应力σn以拉为正,切向应力τs以逆时针为正,而其大小与变形成的关系为(不考虑变形耦合效应):式中,ks、kn分别为切向和法向刚度系数,可分别近似取λ+2G、G,其中λ、G分别为拉梅常数、剪切模量,并按较强一面取值。相对位移可由整体坐标系下不连续界面两侧桩和土节点单元的位移差并经坐标转换后求得:式中:T为坐标转换矩阵;up、vp为整体坐标系下桩域节点x向、y向位移;us、vs为整体坐标系下土域节点x向、y向位移。由式(6)可将桩域或土域中某一节点的位移向量up、us表示如下:式中:ΦpT(x)、ΦsT(x)分别桩域和土域中节点形函数;u为支持域内各节点位移所组成的向量。由式(7)―式(9),采用变分原理可得界面元对总刚度矩阵的贡献为(对荷载列阵贡献为0):将式(10)展开并采用高斯积分整理后叠加至系统整体平衡方程,可近似考虑非变形协调的桩土接触分析,具体分析可参见图4所示的接触算法。3桩身失稳临界荷载的无单元分析模型及求解基于上述处理方法编制了基桩非线性屈曲分析无单元伽辽金法程序(ProgramforBucklingAnalysisofPilesbyEFGM),主要步骤为:1)在桩土分析域Ω内及边界Γ上分布若干节点(即场点,可均匀或不均匀分布,位移变化梯度大或应力集中处可加密节点分布),构成无单元分析模型,并在该模型上给定荷载与边界条件。2)选择合适基函数、权函数及影响域等,并基于上述场点分布,用移动最小二乘法(MLS)获得位移插值函数或形函数,由此构造近似场函数。3)根据系统控制微分方程及边界条件,用步骤2)得到的形函数及系统控制微分方程的强或弱积分形式,建立基于整体坐标系下的离散方程组。4)由系统方程组求解桩身失稳临界荷载。求解时采用基于Crisfield-修正Riks-Ramm弧长法的非线性屈曲分析。为更准确预测桩身失稳极限荷载,初始载荷比例系数不能设太大,最大弧长放大系数也应有所限制,否则后续加载步长太大而对极限荷载分析不利;但也不能太小,避免过大增加总增量步,浪费计算时间。因此,应经试算或按经验合理取值,使得失稳前后增量步弧长足够小而其他加载阶段的弧长则可适当放大。本文采用的办法是:由桩身线性屈曲分析结果估计桩身失稳发生的大致载荷区间,并以此为界将加载过程定义为若干个历程,合理设置各历程中的控制参数。4工程应用4.1桩身结构及试验方案茅草街大桥位于湘北洞庭湖软土地区,主桥为三跨连续自锚中承式钢管混凝土拱桥,采用超长钻孔嵌岩灌注桩基础。为探讨这类桩的承载特性,在主墩附近对一试桩进行了竖向和水平载荷试验。场地内的岩土层分布及其力学指标,如表1所示。试桩混凝土为C30,直径1.0m,长61.5m,入土深60.0m,地面上长1.5m部分为加载所需桩头,桩端进入弱风化泥质砂岩,正反循环旋转钻成孔,泥浆护壁,为防塌孔上部设钢护筒(外径1.3m,壁厚7mm,埋深10m)。成孔后至钢筋笼下放前钻机不停,压浆清理沉渣,混凝土浇注前实测沉渣厚约50mm。桩身内共埋设了19组、38个钢筋应变计和34组、136个混凝土应变计。待桩身混凝土强度达设计要求(桩身下、桩身中、桩身上三部位试块平均抗压强度分别为31.5MPa、27.2MPa和40.6MPa)后,先后进行了竖向和水平载荷试验。其中,竖向载荷试验采用锚桩-反力梁加载,锚桩直径1.8m,反力横梁宽1.7m、高2.7m,两端通过圆立柱与锚桩相连。采用行程60cm的9个250t千斤顶联动加载。4.2桩身荷载变化限于篇幅,这里仅给出桩身轴力和侧阻与桩土相对位移关系曲线(图5和图6)。加载至16200kN时桩顶沉降达32.38mm且能稳定;因此时加载量已超过预估极限承载力,故将原定加载步减半,但加至16720kN时,荷载只能勉强维持且沉降出现增大趋势;继续加荷至17280kN并维持约5min后,沉降增加速率明显减小,并随之从地下发出一声很闷的巨响,试桩随即偏位,地面出现裂缝,沉降无法稳定,千斤顶油压也骤降至11880kN后才勉强稳定。与此同时,离桩底约15.0m范围内的混凝土应变计读数也从有到无。所有这些表明,桩身产生了压屈破坏,失稳极限荷载可取为16720kN。4.3in法数值分析先按规范查表法进行计算。由水平载荷试验结果(Hcr=175kN,Hu=325kN)反算出桩侧地基土当量m值约为4.0MN/m4,由此按m法求出桩土变形系数α=0.3739m-1。然后,由桩身无量纲长度αh及桩端嵌固查表计算桩身稳定计算长度lp=8.54m,进而求得Pcr=133150kN。显然,这一结果远大于实测值,主要原因在于难以从桩土复杂作用体系中取定准确而又简单的参数。为此,采用本文建议的Galerkin法进行数值分析,并先假定桩土接触面为弹性小变形下的变形协调,按第一类稳定问题进行分析。如图7所示,计算域取径向和桩端以下土域范围均为10倍桩径,边界条件取为简支;桩域和土域分别采用375个和5740个“点单元”进行离散,接触面上共布置245个“接触点对”,且土域内的点单元呈内密外疏分布。桩周土体的参数参按表1取值,而桩身混凝土计算参数则按表2取值。然后,进行了线性与非线性增量屈曲分析,求得Pcr分别为22910kN、19300kN,其与实测值相对误差分别为37.0%和15.4%,误差较大的原因在于:弹性分析无法反映实际桩土体系弹塑性破坏结果;其次是计算时没有考虑桩土材料非线性及接触面上的滑移、开裂等非变形协调。因此,继续考虑桩土材料与接触非线性,采用图4非线性迭代算法进行分析。计算时桩土参数分别按表2与表1取值。为简化计算,桩土接触面摩阻系数统一取加权近似值0.5。然后,以非线性屈曲分析获得的19000kN为拟施加最大荷载,并经试算设定初始增量步荷载与总荷载比为0.01,最大增量步荷载与总荷载比为0.5。计算表明,当加荷至17650kN后,桩顶计算沉降量增幅较大(如图8),且反复迭代计算难以收敛,由此可判定桩顶荷载已达一最小稳定极限,并可认为该值即为桩身失稳荷载。尽管桩身失稳荷载计算与实测值误差较小(相对误差5.6%)。但图8所示计算与实测荷载-沉降曲线仍存在一定差别,这可能是由于没有考虑加载偏心及桩身可能初始缺陷等因素。5桩土接触分析(1)基于第二类稳定问题引入新型无单元伽辽金法建立的高承台嵌岩灌注桩屈曲分析数值模型更便于考虑桩土共同工作中的材料非线性