李 磊
(中铁二十二局集团有限公司北京工程勘察设计院，北京 102600)
对于桥梁常见的桩基础，现行的设计方法多采用杆系级地基模型。土体多采用弹簧单元，考虑土的地基系数这样一个参数来确定反力分布，以此来表达土对结构的作用。这是基于E.Winkler地基模型的假定，符合现今的桥梁设计规范以及被学术界普遍认可的研究成果。其计算方法和结果应用广泛，也符合相关的理论和相应的试验结论，在大多数情况下和一定的条件或者假设下，是成立的，广泛应用于实际的工程设计和理论研究[1]。然而对于一些地质情况比较特殊或者比较复杂，或者重要的、大型的桥梁结构，以及进一步的理论研究中，仅仅采用简化的E.Winkler地基模型显然是远远不够的。
建立实体级的土体模型显然对于判断规范计算的合理性和精度或者进一步研究土体与结构的作用是必要的。土的理想化处理的最简单形式是介质土的连续统一化，即设土为连续介质和具有线弹性性质(隐含线弹性和可塑性假设)。随着计算手段的提升，大型的有限元计算软件如ANSYS、SAP2000等的引进和广泛应用，大型的、复杂的非线性分析成为可能，使得土和结构相互作用的非线性分析越来越多地进入到人们的研究视野[2]。本文就是探求桩土接触分析的有限元仿真及ANSYS实现方法。
ANSYS中的土体是采用DP准则建立的一种模型，DP准则是摩尔-库伦准则的近似，是在密塞斯(Von Mises)准则的基础之上又考虑了平均主应力对于土体抗剪强度的影响发展起来的广义密塞斯(Von Mises)准则。结构的屈服面不随材料屈服而发生逐步变化，没有固定的强化发展准则，结构的本构模型可以近似的看作理想的弹塑性体，可以考虑关联流动法则[3]。采用该准则确定的屈服强度考虑了侧限压力，考虑了体积膨胀(由屈服引起)，没有考虑温度影响。广泛适用于土、岩体、混凝土等材料，其主要的影响因素为内摩擦角、剪胀角和粘聚力。这些参数的选择可以参看相关的文献资料、国外的设计规范。在一些深层次的理论研究中，往往也可以采用参数识别或者试算的办法，利用一些勘察资料初步确定土体参数，然后通过计算结果与试验结果修正各项参数，直到二者在一定范围内符合，然后再大规模的分析与计算。
接触分析是高度的非线性行为，在确定结构的受力状态之前，并不能确定接触区域和表面之间是否仍然接触或者突然变化，是否需要接触摩擦[4]。其决定因素往往很复杂，往往不同荷载、不同材料、不同的边界条件、不同的摩擦模型或者其他一些影响因素，都对其有影响，而这些往往都是非线性的。采用ANSYS结构分析软件模拟该部分受力状态时，通常采用点-点、点-面和面-面三种接触方式，这里主要采用面-面接触来研究桩土接触分析。面-面接触通常将刚性面作为目标面，柔性面作为接触面。前者可以模拟桩接触面，后者可以模拟土接触面，这样定义接触对，来模拟桩土接触部分的力学行为，从而分析桩身或者土体(接近桩基础部分)的受力和变形，达到分析桩基础的沉降和桩土应力分布的求解目的。
ANSYS的接触分析过程主要包括：建立模型、划分网络、识别接触对、定义刚性目标面和柔性接触面、设置实常数与关键字、控制刚性目标的运动、设置边界条件、定义求解分析选项、定义荷载步选项、求解并检查结果几个步骤。
3.1 问题的描述
设计一个简单算例，分析某混凝土方桩在有水平作用情况下桩、土的应力分布。通过求解结果，将数据进行一定处理，并与桩采用梁单元、土体采用E.Winkler地基模型的弹簧单元来模拟求解出的桩身内力、桩基沉降的结果进行比较。如图1所示方桩长10 m，入土深度8 m，截面边长2 m，在桩基础的顶部，荷载作用为：弯矩M=1 200 kN·m，水平力H=150 kN，轴力N=2 100 kN。
图1 方桩计算图示
3.2 计算参数的选取
桩身混凝土弹性模量Ec=2.6×107kN/m2，泊松比μc取0.167，密度ρc=2 500 kg/m3。土体的弹性模量为E=0.95×104kN/m2，泊松比μ取0.42，密度为ρ1 900 kg/m3。土体粘聚力为19 kPa，摩擦角为24°，膨胀角取20°，土与桩的摩擦系数为0.3。
采用m法时，m值取2 000～15 000 kN/m4试算，根据土的性质分别取3 000 kN/m4、5 000 kN/m4与接触分析结果进行对比。
3.3 建模思路
将受力结构的桩基础和周围土体均采用SOLID45单元模拟，桩身按实际尺寸建模，土体在深度方向上取桩的3倍，在x～z平面方向取5～8倍。由弹性模量的量级，指定桩为刚性体，土体为柔性体。将桩基础上的接触面假定为刚性目标面，用TARGE170(3D)单元，土体用CONTA173或CONTA174模拟，假定为柔性接触面。在本工程实例中，采用CONTA173(4节点)模拟。
3.4 计算结果与比较
如图2～图4分别是ANSYS计算求解出的桩的应力分布和土的位移、应力分布云图，这与预期结果和试验或者理论结果是相符的。采用与前面的算例相同的计算参数，对ANSYS的每个截面提取应力后，求解出桩身内力和桩顶沉降，与m法的计算结果进行比较，如图5～图7所示，m值取值3 000 kN/m4或5 000 kN/m4，其内力计算结果基本相同。
图2 桩体计算结果 图3 土体表面视图下x向应力
图4 土体表面视图下z方向计算结果
图5 不同方法桩身位移比较
图6 不同方法桩身剪力比较
图7 不同方法桩身弯矩比较
由结果可以看出，m法的计算结果与桩土接触分析的计算结果基本相同，差异主要由土参数选取引起，桩和土的应力分布与单桩的分布情况类似，预期结果与试验或者理论结果是相符的。
4.1 问题的描述
群桩结构和各项参数如图8所示，其中荷载N=11 291.69 kN，H=836.23 kN，M=7 782.9 kN·m。自由长度l0=12.81 m，入土深度h=25.19 m，最小桩中心距l中=3.3 m，混凝土重度为15.0 kN/m3(已扣除浮力)。
图8 桩基础结构
4.2 计算参数的选取
设计桩身混凝土材料的弹性模量Ec=2.8×107kN/m2，泊松比取0.3，密度ρc=2 500 kg/m3。土体弹性模量为0.95×104kN/m2，泊松比μ=0.2，密度ρ=1 900 kg/m3。土体粘聚力为19 kPa，摩擦角为24°，膨胀角取20°，土与桩的摩擦系数为0.3。
4.3 建模
群桩的建模和单桩相同，将受力结构的桩基础和周围土体均采用SOLID45单元模拟，桩身按实际尺寸建模，土体在深度方向上取桩的3倍，在x～z平面方向取5～8倍。由弹性模量的量级，指定桩为刚性体，土体为柔性体。将桩基础上的接触面假定为刚性目标面、用TARGE170(3D)单元，土体用CONTA173或CONTA174模拟，假定为柔性接触面。在本工程实例中，采用CONTA173(4节点)模拟。
4.4 计算结果与比较
不同算法的计算结果如表1所示。
表1 算例中的结果比较
由结果可以看出，m法的计算结果与桩土接触分析的内力计算结果基本相同。桩身位移较小，这是由于实体模型的参数选取产生的差异。桩土接触分析结果的沉降比m法计算结果稍大，但是同样因为与规范的假设不同，结果是不同的，这与预期结果和试验或者理论结果是相符的。采用与前面的算例相同的计算参数，对各项计算结果比较，对ANSYS的每个截面提取应力后，求解出桩身内力和桩顶位移，与m法的计算结果进行比较，其结果具有一定的参考价值。图9和图10分别是桩身的弯矩和剪力分布图。由桩土接触分析的结果可以看出，只要正确的确定了土体的参数，也能得到较为精确的结果。
图9 桩身弯矩比较 图10 桩身剪力比较
本文通过对桩-土接触分析的研究探讨，结合典型单桩、群桩基础的算例，并将桩土接触分析结果与弹性地基梁法结果进行比较，得出结论：只要DP模型的参数选择合理，两者结果一致，特别是桩身内力分布、桩身横向位移等一些工程中最需要得到的结果，m法和ANSYS方法的计算结果可靠，精度也很高；
但是对于基础沉降、地基以下较深处土的反力分布等有待进一步的研究。