接触问题广泛存在于土木、水利等各类工程问题中．在桩底位置，桩与土的相互作用是一种典型的法向接触问题．随着计算机技术的快速发展，以有限元方法[1]为理论基础的各种商业软件程序(ANSYS，Abaqus，ADINA，LS-DYNA等)被开发应用．人们可以对工程中各种结构问题的应力场和应变场进行定性和定量的计算[2-3]．目前，商业软件的接触算法主要包括罚函数法[4]、拉格朗日乘子法[5]、增广拉格朗日乘子法[6-8]等．然而，这些算法在处理接触问题中仍然存在一些不足．罚函数法实现简单，不增加新的未知数，是一种工程近似方法．在实际工程计算中，通过引入穿透值(接触位移)和接触刚度，允许接触体之间出现微小嵌入，进而确定接触力的大小．一方面，选取的罚参数对接触力的收敛性和精度有很大影响，接触刚度设置过小会导致虚假穿透过多，过大又会导致刚度矩阵的病态性；另一方面，对桩采用实体单元进行装配时，除了增加额外的建模工作量以外[9]，桩土体边缘单元划分的近似性与装配过程的不符合性还会导致两种结构间发生互相排挤，从而在初始应力场形成较大的虚假应力．拉格朗日乘子法能精确满足接触条件，但是引入多余自由度会增加系统方程的维数和求解难度；接触状态突变会导致接触力的突变，产生接触状态的振动式交替改变；还可能出现对角线零主元情况，对计算求解器的要求较高．Simo等[6]提出了摄动拉格朗日乘子法；Simo等[8]、王水林等[7]提出了增广拉格朗日乘子法．增广拉格朗日乘子法求解精度同样依赖于罚参数值的大小，与罚函数相比减少了刚度矩阵的病态性，增加了迭代次数；与拉格朗日乘子法相比，没有对角线零主元情况，但其穿透值不为0．互补问题的最显著特征是含有互补性条件，即要求两组非负变量(或变量的函数)对应分量的乘积为零[10]．Costa等[11]提出了一种非单调混合互补问题来解决单边摩擦接触稳定性问题．Kwak等[12-13]将三维摩擦问题转化为非线性互补模型，推导了相应的非线性方程组，并用同伦算法求解．Bathe等[14]将互补函数用约束函数进行逼近来处理接触问题，提高了收敛性．孙苏明等[15]提出摩擦接触弹塑性分析的数学规划法．陈国庆等[16]建立了三维弹性摩擦接触问题的互补变分原理及有限维非线性互补模型，并提出极小化算法．陈万吉等[17]推导出了接触问题的互补虚功原理，互补能量原理及变分不等式等变分提法，并讨论了三类提法的关系．郑宏等[18-19]将DDA的接触问题归结为一个非线性混合互补问题．林珊等[20]就面-面接触模型通过调节积分点数目对求解精度进行控制，并在此基础上建立二维接触有限元模型．本研究提出一种混合结构单元和实体单元的端承桩法向接触互补问题的非光滑牛顿算法．1 法向接触的非线性互补模型实际工程中主要在端承桩底部位置考虑法向接触，如图1所示．图中：g为端承桩底部的桩土间隙；nos为桩土底部接触位置局部坐标．10.13245/j.hust.239216.F001图1桩底接触模型法向应力σn=F/A，桩土间隙g=up-us，式中：A为端承桩的横截面积；F为轴向接触力，以压为负；up为桩体的竖向位移；us为土体的竖向位移．σn和g满足Karush-Kuhn-Tucker(KKT)条件：g≥0;-σn≥0;gσn=0, (1)利用非线性互补函数可以将式(1)转化为等式问题，从而把端承桩的法向接触问题归结为解方程组的问题．选择不同的互补函数，求解效果不同．∀ x，y∈R，Fischer[21]使用如下函数形式，即Fischer-Burmeister函数ψ(x,y)=x2+y2-x-y．(2)Fischer-Burmeister函数只在原点不可微，具有很强的半光滑性，经常用来转化互补问题[20]．因此，可以利用式(2)将式(1)重新表述为非光滑方程σn2+g2+σn-g=0, (3)式(3)称为等式约束平衡方程．2 法向接触的非线性有限元方程组在图1中，桩体采用梁单元，土体采用实体单元，根据桩土间的受力情况，端承桩与土体之间的有限元平衡方程如下∫VBTσdV-Fext-σnANTn=0, (4)式中：B为梁单元与实体单元装配后的总的应变矩阵；σ为模型系统内对应高斯点应力；V为系统总体积；Fext为系统外力；N为总体形函数矩阵；n为法向向量．将式(3)与(4)联立，可到端承桩法向接触的非线性有限元方程组∫VBTσdV-Fext-σnANTn=0;σn2+g2+σn-g=0. (5)将式(5)的2个算式分别记为f1(U，σn)和f2(U，σn)，U为总体位移矩阵．g的离散表达式为g=nTNU．令x=(U，σn)T，f =(f1，f2)T，式(5)可写成fx=0．(6)式(6)可采用牛顿-拉弗森算法[22]进行求解，其迭代序列如下xk+1=xk-f '(xk)-1f(xk)，式中雅可比矩阵f ' (x)可以写作f '(x)=∂f1∂U∂f1∂σn∂f2∂U∂f2∂σn=K11K12K21K22，其中，K11为梁单元与实体单元装配后的刚度矩阵；K12=-ANTn；K21=gg2+σn2-1nTN；K22=gg2+σn2+1．本文算法采用的收敛准则x(k+1)-x(k)≤ε，式中ε为容许误差．3 算例分析3.1　单根端承桩算例将长度100 m的单根端承桩埋入高300 m、宽200 m的土中，埋深70 m，并位于土体中央．桩的泊松比为0.2，弹性模量为21 GPa，圆形横截面直径为1 m，抗弯惯性矩为0.25π m4；土的弹性模量为20 MPa，泊松比为0.3，重力加速度取9.8 m/s2，土体密度为2 600 kg/m3，黏聚力为10 kPa，摩擦角为22°．在桩体埋入土体达到稳定后，于桩体顶部受到0.05 m(方向向下)的初始位移，土体两侧采用法向约束，土体底部固定．由于是端承桩，因此不考虑桩土间的侧向摩阻力，对桩侧位置桩土之间施加了横向位移协调条件，且仅由上部荷载引起位移增量和应力增量．桩体划分为50个平面梁单元，土体划分为2 400个四节点四边形单元．对桩体分别采用实体单元和结构单元，通过商业软件Abaqus直接求解，其中实体单元桩底与土体之间采用接触单元，结构单元与土体在桩底固结在一起．选取桩底处水平向右位置的位移与应力数据组，将计算结果进行对比．土体沿桩底位置(桩底中心水平距离D)的水平位移(L)对比如图2所示，可以看出：本文算法计算结果和另外两种方法拟合度较好；由于单元划分的关系，采用实体单元在近桩底附近位置计算结果的平滑性较差．本文结果不仅更接近采用实体单元的计算结果，而且更加光滑．10.13245/j.hust.239216.F002图2土体沿桩底位置的水平位移对比土体沿桩底位置的竖直位移(H)对比如图3所示，采用桩体实体单元的计算结果是比较准确的，本文算法的计算结果与其非常接近，由于不考虑接触面单元，因此在计算效率和收敛性上具有较大优势；商业软件采用结构单元计算时，由于看作是简单的点-面的接触，对周围土体的影响能力被削弱，因此计算结果偏小．10.13245/j.hust.239216.F003图3土体沿桩底位置的竖直位移对比在近桩底位置附近3种方法的总位移差别不大，在接近土体边界附近采用Abaqus的结构单元计算得到的结果偏小，这可能是由于其点-面接触算法对远端土体的位移扰动较实际情况偏小导致的．竖向应力(p)的计算结果对比如图4所示．各算法得出的近桩底位置的竖向应力差别较大．这是由于商业软件采用桩体实体单元计算时在桩底位置发生应力集中，竖向应力产生过大；采用桩单元计算时由于采用局部Embedded的接触方式影响范围过小，导致桩底产生的竖向应力过小；本文算法的计算结果介于商业软件两种单元方式之间．在远离桩底位置，三者的计算结果都比较接近．10.13245/j.hust.239216.F004图4竖向应力的计算结果对比3.2　端承桩群桩将长度100 m的11根端承桩埋入高300 m、宽200 m的土中，埋深70 m，并位于土体中央．埋入端承桩的泊松比为0.42，弹性模量为21 GPa，圆形横截面直径为1 m，抗弯惯性矩为0.25π m4，上部横梁弹性模量远大于埋入桩，为210 GPa．土的弹性模量为20 MPa，泊松比为0.3，重力加速度取9.8 m/s2，土体密度为2 600 kg/m3，黏聚力为10 kPa，摩擦角为22°．在桩体埋入土体达到稳定后，于桩体顶部受到1 m(方向向下)的初始位移，土体两侧采用法向约束，土体底部固定．由于是端承桩，因此不考虑桩土间侧向摩阻力，对桩侧位置桩土之间施加了横向位移协调条件，且仅由上部荷载引起位移增量和应力增量．桩体划分为600个平面梁单元，土体划分为7 200个四节点四边形单元．土体沿桩底位置的水平位移对比如图5所示，两者水平位移先增大后减小，趋势类似．近桩底到土体中部位置，本文算法得到曲线较为平缓，光滑性好，对于土体变形的横向传递较为全面；远端位置两者的符合度较高，规律性较为一致．10.13245/j.hust.239216.F005图5土体沿桩底位置的水平位移对比土体沿桩底位置的竖直位移对比如图6所示，商业软件与本文算法的计算结果相近．总位移的整体趋势都是远离桩底中心位置逐渐减小，最终趋近于0 m．10.13245/j.hust.239216.F006图6土体沿桩底位置的竖直位移对比竖向应力与商业软件计算结果对比如图7所示，整体拟合情况较好，在群桩底部位置区域(0~50 m范围)，商业软件采用实体单元计算的应力集中现象较为明显，结果偏大．在较远位置，两者的计算结果都十分接近．10.13245/j.hust.239216.F007图7竖向应力与商业计算结果对比可见，本文方法的计算结果与经典算法基本一致，但更光滑且计算量较小．4 结语针对现有商业软件在求解端承桩法向接触问题中的不足，提出了相应的混合结构和实体单元的非光滑牛顿算法．给出了单桩和群桩的算例，并将计算结果与商业软件Abaqus的计算结果进行了比较，结果表明本文算法的计算结果与经典算法基本一致，但更光滑，说明所提算法是可靠的．