Boulanger桩-土非线性动力相互作用p-y模型二次开发的验证

渠继鑫，周理安

（1.天津城建大学 土木工程学院，天津 300384；2.北京建筑大学土木与交通工程学院，北京 100044）

桩-土相互作用（pile-soil-interaction，简称 PSI）属于桩基抗震设计中一个重要的课题，刘立平等[1]数值分析了桩-土-结构相互作用对结构弹塑性变形特性的影响，结果表明桩-土相互作用对结构的弹塑性变形的影响不容忽视.在桩土动力相互作用分析中，Winkler地基梁模型得到了普遍应用.这种方法将桩周土对桩的动力阻抗用连续分布的相互独立的弹簧和阻尼器替代，具有简便实用，物理概念清楚，计算量小等优点.具有代表性的模型包括Matlock模型[2]、Kagawa 模型[3]、Novak 模型[4]和 Nogami模型[5].但以上模型存在模型粗糙，不能考虑非线性，不能反映桩土界面分离等缺点.孔德森[6]基于Nogami模型，通过增加等效模型的自由度和相应参数，改善了实际桩基动力阻抗的逼近效果.姜忻良等[7]进行了群桩-土-偏心结构相互作用体系的振动台实验；凌贤长等[8]采用振动台实验数据，对建立可液化场地桩-土-桥梁结构地震相互作用数值模拟的二维分析模型的有效性进行验证分析；Yang[9]建立了较为理想的自由场土侧本构模型，考虑自由场土对桩土相互作用的影响.Juirnarongrit[10]通过在桩头处添加旋转约束弹簧来模拟桩帽约束对桩土相互作用的影响，并对p-y动力法加以分析和评估.Malhotra[11]系统地介绍了地震下桩-土相互作用的关键因素，并对每个因素进行了分析.Brandenberg[12]进行了液化土壤中群桩下桩-土相互作用的离心实验.此外，刘立平等针对常用的刚结模型、铰结模型、接触模型、建华模型和弹簧模型5种桩-土动力相互作用分析模型，对模型模拟的有效性进行对比验证分析，结果表明，不同的分析模型有不同的适用范围.

本文基于Boulanger提出的p-y模型，通过对现有计算软件的二次开发，发展了桩-土相互作用的有限元计算程序，并采用单向推覆和低周反复土体动应力与动应变实验对程序的有效性进行验证.

1　p-y模型

Boulanger于1999年提出了一种模拟[13]桩-土相互作用的p-y模型.该模型由弹性单元（p-ye）、刚塑性单元（p-yp）、开合缝单元（p-yg）和黏滞阻尼单元4部分组成.LS-DYNA提供了完善的自定义材料接口，使用户方便进行自定义本构模型的开发.本文利用Fortran语言将p-y非线性本构模型编入主程序中.编程需要用到LS-DYNA主程序的接口变量有：应变增量eps（）、上一时间步应力sig（）、上一时间步历史变量hsv（）、当前时间步dt1和材料参数阵列cm（）.

2　模型验证

2.1　实验概况

采用GCTS空心柱扭剪仪开展了一系列不排水条件的动三轴循环荷载实验，验证不同的固结围压、动应力幅值和加载波形等工况下二次开发程序的有效性.利用单向推覆实验和动三轴循环荷载实验进行验证需要对二次开发程序进行修改，由于该实验不能模拟桩土分离及桩侧土体的侧向摩擦阻力的效果，数值模拟应剔除开启弹簧产生的横向摩擦阻力的作用，考虑单侧土体弹性弹簧、刚塑性弹簧、关闭弹簧的作用，土极限力系数的1.8应去掉.仪器在循环荷载、外围压的控制下，使土体单元的受力特征更加接近土体的实际受力特征，循环荷载模拟了桩对土体的作用，围压模拟了深度对土体性质的影响.正向幅值由于桩挤压土体的作用而产生，模拟桩单向作用下二次开发模型弹性和刚塑性弹簧的性能；负向幅值由于围压对土体的作用而产生，模拟了桩开始反向移动尚未与土体分离时的相互作用，模拟桩单向作用下二次开发模型刚塑性弹簧的性能；由负向幅值向正向幅值的过程由于桩再次正向移动产生，模拟桩二次开发模型关闭弹簧的性能.实验设计取值见表1，土体基本物理性质[13]见表2.

表1　实验方案

表2　土的物理性质

2.2　单向推覆实验

图1为固结围压75 kPa和150 kPa时静三轴实验结果和数值模拟结果对比曲线.可以看出土体的极限应力和初始弹性模量随固结围压的增大而增大.从图1中可以看出，所开发的p-y模型在土体处于弹性阶段时能精确地模拟桩-土相互作用过程，随着应变的增加，固结围压为75 kPa的工况模拟精度较固结围压150 kPa工况的精度有所降低，原因是低围压下土体的静三轴实验具有较大的离散性，实验过程较难控制，表现为数值模拟土体较实验土体的刚度退化慢，而高围压下数值模拟土体和实验土体的刚度退化吻合更好，数值模拟结果随应变的增加不断逼近于极限应力.总体而言，所开发的p-y模型在土体单向推覆分析中具有很高的计算精度.

图1　固结围压75 kPa，150 kPa静三轴对比曲线

2.3　循环荷载实验

土的动应力-应变关系是土体动力分析时必不可少的基本土性关系.土体在动荷载作用下的实验表明，在周期荷载下，土的动应力-应变关系具有3个基本特点：非线性、滞后性和应变累积性.以下将针对前述实验结果从动应力幅值、固结围压、加载波形等方面对模型的有效性进行分析.图2-4为固结围压75 kPa，荷载幅值分别为 20，30，40 kPa时动三轴实验结果和数值模拟结果对比曲线；图5-7为固结围压150 kPa，荷载幅值分别为 20，30，40 kPa时动三轴实验结果和数值模拟结果对比曲线.

图2　固结围压75 kPa，荷载幅值20 kPa动三轴对比曲线

图3　固结围压75 kPa，荷载幅值30 kPa动三轴对比曲线

图4　固结围压75 kPa，荷载幅值40 kPa动三轴对比曲线

2.3.1动应力幅值的影响

由图2-7可以看出，土体的动应变随着振动次数的增加而不断增长.当动应力幅值较小时，围压75kPa和围压150 kPa的试样均不会达到破坏状态，能够振满1 000圈停止，加载初期土体动应变增长较快，振动后期动应变增加的速率变慢，土体被逐渐压密.当动应力幅值较大时，试样加载初期土体力学性能稳定，动应变随振次的增加变化较小，但达到一定振次时，试样动应变增长速度变快，土体的刚度出现一定的退化，且滞回环面积增大.动应力幅值为20，30，40 kPa工况的数值模拟结果均能较好地模拟土体的应变累积和滞回耗能，同时也能较好地模拟土体刚度退化到初始刚度80%时的破坏.由于高频次振动的实验过程中仪器对应力幅值的控制有所偏差，而数值模拟按照精确的预设应力幅值加载，所以各工况数值模拟结果的应力幅值和滞回面积均高于实验结果的应力幅值.

2.3.2固结围压的影响

由图2-7可以看出，在同样振次下，高围压下土体动应变增长速度较低围压下土体动应变增长速度快，且在应力幅值较大时，高围压下土体达到破坏时的振次较低围压下土体达到破坏时小.同时，高围压下土体刚度较低围压下土体刚度大，相同应力幅值下滞回环面积小.固结围压为75 kPa和150 kPa工况下，极限应力的增长、土体刚度的提高以及高围压下滞回面积的减小等，数值模拟结果与实验数据都吻合较好.就比较而言，固结围压为75 kPa的工况数值模拟结果精度较固结围压150 kPa工况偏低，原因是低围压下数值模拟土体较实验土体的刚度退化慢，使模拟的低围压土体刚度普遍高于实验值.

图5　固结围压150 kPa，荷载幅值20 kPa动三轴对比曲线

图6　固结围压150 kPa，荷载幅值30 kPa动三轴对比曲线

图7　固结围压150 kPa，荷载幅值40 kPa动三轴对比曲线

2.3.3加载波形的影响

由图2-7可以看出，方波加载时的动应变较正弦波加载时的动应变大.这是因为方波的施力状态是正向最大加载，之后立即反向施加最大负载，使土体承受的剪应力方向改变.数值模拟结果很好地反映出土体在方波工况较正弦波工况耗能大、破坏振次少、刚度退化快等特点.方波工况数值模拟结果精度较正弦波工况的有所降低，原因是方波能量较大，相应土体产生的不可逆塑性应变和刚度退化均较大，土体产生的非线性增强.

3　p-y模型算例

利用LS-DYNA模拟一个桩基础桥墩作为算例.模型墩高23 m，桩长34 m并划分为34个单元，相应的有34个p-y单元.土层按土质划分为6层，从上到下依次为淤泥、淤泥质土、粉质黏土1，粉质黏土2、粗砂、粉质黏土3.p-y单元一侧与桩相连，在p-y单元另一侧输入各自由场土层加速度时程.模型示意图如图8所示.分别提取了2 m淤泥土土层、10 m粉质黏土土层，桩侧p-y单元的应力应变曲线，如图9-10所示.

图8　桩基础桥墩模型

图9　深度2 m下p-y单元应力应变曲线

图10　深度10 m下p-y单元应力应变曲线

4　结论

本文以LS-DYNA二次开发为基础，结合该软件的用户可编程特性（UPFs），编写并导入了考虑桩-土相互作用的p-y动力非线性本构模型，对模型的有效性进行了实验验证.主要得出如下结论：

（1）通过土体静三轴实验及动三轴循环荷载实验表明，所开发的模型具有较高的模拟精度，能同时模拟不同土质下桩-土非线性动力相互作用关系.

（2）非线性p-y模型在土体基本物理性质的基础上，考虑了桩径、埋深、黏滞阻尼、土体弹性及刚塑性区域、土体回弹、桩土分离、闭合等因素，能够较全面地模拟桩侧与土体之间的相互作用.

（3）所开发的模型简化了桩-土相互作用数值建模的复杂性，减小了桩-土接触面处理和软土特性参数选取带来的误差，并且大大节省了运算速度.