弹性地基Timoshenko梁单元在ABAQUS软件中的应用

杨钊，许建聪，余俊

（1.同济大学岩土及地下工程教育部重点实验室，上海 200092；2.中南大学土木建筑学院，湖南 长沙 410075）

弹性地基Timoshenko梁单元在ABAQUS软件中的应用

杨钊1，许建聪1，余俊2

（1.同济大学岩土及地下工程教育部重点实验室，上海 200092；2.中南大学土木建筑学院，湖南 长沙 410075）

利用ABAQUS软件中的自定义单元接口，采用Fortran语言开发弹性地基Timoshenko梁单元程序.通过与经典解的比较，验证所编弹性地基Timoshenko梁单元程序的准确性.该单元不仅考虑了地基弹簧受拉脱开的特点，而且还考虑了曲形梁单元内结点不在一条直线上的特点.采用所编写弹性地基梁单元分析青草沙原水过江输水工程，计算规律与已建类似工程实测结果相同.

弹性地基；Timoshenko梁；ABAQUS软件；单元子程序

在地下结构的计算领域，已有许多用于考虑结构-土体相互作用的计算方法，但是荷载结构法仍是目前使用最广泛的一种方法.为了在研究周围地层对结构的约束作用及地层对地下结构的反作用力时考虑地下结构变形，将Winker弹性地基理论引入到荷载结构法中［1］.同时，提出弹性地基上的Timoshenko梁，以考虑横向剪切变形对厚梁的影响.目前，大量商品化软件如MARC，SAP等均包含弹性地基梁单元，但均没有考虑地基弹簧在受拉时的脱开情况.ANSYS，ABAQUS等大型通用有限元法计算分析软件中未包括弹性地基梁单元，只能通过在梁单元结点上加弹簧单元来近似模拟弹性地基的作用［2-3］.这种做法只能在梁单元足够小的情况下，才能近似等价于弹性地基作用［4］.本文以ABAQUS软件为平台，开发基于Winker地基理论与Timoshenko梁理论的地基梁单元.

1 Timoshenko梁单元刚度阵与荷载阵

设一Timoshenko梁的弯曲刚度为EI，剪切刚度为S，压缩刚度为EA，长度为L，梁上有侧向分布荷载q（x），径向集中力Pi，集中力偶Mj，轴向集中力Fk.因此，其总势能［5］为

Timoshenko梁单元的基本特点是挠度ω、轴向压缩量u和截面转动θ各自独立插值.即

式（2）中：n为单元的结点数；Ni是Lagrange插值多项式.将式（2）代入式（1）中，由δΠ=0可以得到有限元的求解方程为

2 弹性地基梁

弹性地基上，梁在荷载作用下产生变形的同时，地基土也产生了变形.根据弹性地基的局部变形理论，地基土对地基梁的反力集度Pd与地基梁的挠度ω间的关系为

式（4）中：kd为地基反力系数.为了考虑地基弹簧受拉脱开，假定梁的挠度值为负值，kd=0.考虑地基土的应变能后，弹性地基梁的总势能比梁的总势能多出一项［6］，即

式（5）中：b为梁截面的宽度.由式（2）可得

式（6）中：Bd=［Bd，1，…，Bd，n］；Bd，i=［Ni，0，0］.将式（6）代入式（5），取极值即可得单元刚度阵中地基刚度附加项Kde为

梁单元内力计算方程为

式（9）中：N为轴力；Q为剪力；M为弯矩.

3 整体坐标系下的刚度阵

假定单元的方向在整体坐标系下的方向角为β.以单元内任意一点i为例，有

式（11）中：Bi，gol=［Bi，gol，…，Bn，gol］，Bi，gol=Bλii.结合式（6），（7）可得

式（12）中：Bd，gol=［Bd，1，gol，…，Bd，n，gol］，Bd，i，gol=Bd，λii.整体坐标系下，弹性地基梁单元刚度矩阵为

当由前处理所得梁单元的内部结点位于梁单元两端结点的连线上，且内部结点为梁单元的均分点时，由局部坐标与整体坐标的转换关系可得

在实际工程中，特别是以曲梁形式存在的盾构隧道管片，由前处理软件所得梁单元内部结点位于梁两端结点连线之外.针对这种情况，式（14）应为

式（15）中：k为单元在整体坐标系下的斜率.当k=∞时，有

结合式（13）～（16），可以得到在整体坐标系下单元刚度矩阵的表达式.

在数值积分中，若采用精确积分计算剪切变形能项，当梁很薄的情况下，约束条件γ=（dω/dx）-θ不可能在梁单元上处处满足.采用精确积分计算剪切变形能项，将过分夸大剪切应变能的量级而造成剪切锁死.为了避免剪切锁死，对剪切变形项采用缩减积分计算［7］.

4 程序编制

ABAQUS提供用户单元接口子程序UEL，用户通过自定义UEL接口与求解器Standard的接口实现数据传递.UEL有其固定的书写格式与规范，与主程序共享的变量必须在子程序开关予以定义，而主程序通过ABAQUS输入文件（.inp）中的关键字“element，type=Un”来判断是否使用自定义单元.

依据ABAQUS软件二次开发的约定，用户单元子程序UEL应至少包括5部分［8］，分别为：ABAQUS约定的子程序题名说明、ABAQUS定义的参数声明表、用户自定义的局部变量声明表、用户编写的程序代码段和子程序返回与结束语句等.在UEL中，用户需要给出单元的结点数目、结点自由度、材料参数，通过主程序传送给UEL的结点位移及结点位移增量更新单元应力，并最终将单元刚度矩阵（AMA TRX）及单元不平衡力矩阵（RHS）提供给ABAQUS主程序进行迭代求解.

4.1 程序流程

ABAQUS主程序进行迭代求解有如下8个步骤.

（1）计算坐标转换矩阵λi.

（2）计算高斯积分点数与高斯积分点坐标.对于n结点Timensheno梁、轴力项、弯矩项刚度，可由n-1个高斯积分点精确求解；剪切项刚度需要n个高斯积分点才能精确求解.但为了避免剪切锁死的发生，剪切项刚度在计算中也采用n-1个高斯积分点积分.

（3）计算形函数矩阵，并由式（14）～（16）计算积分系数.

（4）由式（11）计算整体坐标系下梁单元的单元刚度矩阵.

（5）由式（9）计算积分点处的轴力、剪力、弯矩.

（6）将主程序传入的位移值代入式（2），（10），计算地基梁单元在积分点处的挠度值.如计算挠度值为负值，则取此积分点处的地基刚度为零，再由式（12）计算地基对梁单元的附加刚度矩阵.

（7）由式（13）计算地基梁单元的整体刚度矩阵.

（8）计算单元的残余力，并判断收敛.如果不收敛，返回ABAQUS主程序进行第i+1次迭代.主程序将根据UEL子程序第i次迭代所得到的单元刚度矩阵与残余力项，计算位移增量与总位移量，然后跳到第（1）步进行UEL的第i+1次迭代计算.

4.2 程序验证

梁荷载和弹性地基梁竖向位移图，分别如图1，2所示.两端自由的弹性地基梁参数：长度l为10m，宽度b为0.5m，高度h为0.5m，梁身的弹性模量为10.0GPa，剪切模量为5.0GPa，地基的刚度系数K为4.0GN·m-3.求梁截面A，B和C的弯矩与挠度［9］.

采用自编弹性地基梁单元计算，将梁划分为20个3结点地基梁单元.计算过程中，可得到每个单元内积分点处的弯矩、剪力、轴力.将积分点处的内力值外推，可以得到单元结点处的内力值.

弯矩和挠度的理论解与数值解的对比，如表1所示.表1中：ωth，Mth分别为弯矩，挠度的理论解；ωc，Mc分别为弯矩，挠度的数值解；σ为相对误差.由表1可知，理论解与数值解相差较小，该误差可能源于有限元数值计算中网格的划分、迭代收敛判断准则，以及其他诸多综合因素的影响.由此可见，编制的有限元程序是可靠的.

图1 梁荷载示意图Fig.1 Loads on the beam

图2 弹性地基梁竖向位移图Fig.2 Vertical displacement of the elastic foundation beam

表1 弯矩和挠度的理论解与数值解的对比Tab.1 Comparison of bending moment and deflection between analytical solution and numerical solution

5 工程实例计算

青草沙原水过江隧洞工程位于上海长江隧道下游约80m处，浦东侧越江点在五号沟，长兴岛越江点在该岛新开河附近，全长7.23km.越江输水管道采用全断面隧道掘进机（TBM）施工，有压输水，设计为圆形断面，衬砌结构外直径为6.8m，管片厚为480mm.考虑冲刷后上部垂直水土压力为439.1kPa，上部水平土压力为313.0kPa，下部水平土压力为375.0kPa，隧道内水压力为404.1kPa.隧道周围地层的地基刚度系数为10.0MN·m-3.采用3根弹簧分别模拟接头的抗弯、抗压与抗剪性能，其刚度系数分别为500（MN·m）·rad-1，5.0TN·m-1，0.5TN·m-1.

计算采用三结点Timeshenko地基梁.为简化计算，在建立地基梁模型的同时也建立三结点Timeshenko梁模型.此梁模型与地基梁模型共结点且结点编号一致，但单元编号不同.将地基梁上的荷载施加到共结点的梁上，取共结点梁单元弹性模量为一极小数.此时，梁单元的存在将简化荷载的施加且对计算结果无影响.由于梁单元结点与地基梁单元共结点，因此，ABAQUS后处理中梁单元的位移场即为地基梁单元的位移场.

图3 衬砌变位矢量图Fig.3 Displacement vector diagram of lining

衬砌变位矢量图，如图3所示.在外周水土压力与内水压力的联合作用下，衬砌结构竖直方向内缩，最大压缩量为2.4mm；水平方向伸长，最大伸长量为2.2mm.衬砌结构的形状由原先的圆形变成扁平的椭圆形.

衬砌截面的弯矩、轴力和剪力图，如图4～6所示.从图4～6可知，管片弯矩的峰值出现在管顶、管底和两腰，其管顶、管底为正值，两腰弯矩为负值.最大正弯矩位于管顶处，其值为102.7kN·m；最大负弯矩位于管腰处，其值为-95.3kN·m.设轴力以受拉为正，受压为负.管片大部分截面受压，轴压值管顶、管底小，而管腰大.轴压最大值位于管腰处，其最大轴压值为150.8kN；管顶部分截面受拉，最大轴拉值为-34.6kN.设剪力以截面呈顺时针转动为正，反之为负，剪力最大值为68.3kN，剪力最小值为-59.7kN.

图4 衬砌截面弯矩图Fig.4 Moment of lining section

图5 衬砌截面轴力图Fig.5 Axial force of lining section

图6 衬砌截面剪力图Fig.6 Shear force of lining section

6 结束语

针对绝大多数商品化软件不能考虑弹性地基梁在受拉时地基弹簧脱开的不足，基于ABAQUS软件平台开发了弹性地基Timeshenko梁单元.该单元不仅考虑了地基弹簧受拉脱开的特点，而且还考虑了曲形梁单元内结点不在一条直线上的特点.采用所编写弹性地基梁单元分析青草沙原水过江输水工程，计算规律与已建类似工程实测结果相同.研究结果表明，所编写的弹性地基梁单元精度高，可供实际工程计算应用.

［1］孙钧.地下结构［M］.北京：科技出版社，1987.

［2］黄群贤，林建华.液化侧扩地基中桩基的有限元分析［J］.华侨大学学报：自然科学版，2004，25（3）：328-330.

［3］贾瑞华，阳军生，马涛，等.既有管线下盾构施工地层沉降监测和位移加载数值分析［J］.岩土工程学报，2009，31（3）：425-430.

［4］杨钊，潘晓明，余俊.盾构输水隧洞复合衬砌计算模型研究［C］//2009全国土木工程博士生学术论坛优秀论文集.长沙：中南大学出版社，2009.

［5］王勖成.有限单元法［M］.北京：清华大学出版社，2006.

［6］朱伯芳.有限单元法原理与应用［M］.北京：中国水利水电出版社，2004.

［7］冯紫良.杆系结构的计算机分析［M］.上海：同济大学出版社，1991.

［8］叶志才，徐磊，王超.基于ABAQUS的三维锚杆单元的开发［J］.三峡大学学报，2008，30（5）：29-32.

［9］黄义，何芳杜.弹性地基上的梁板壳［M］.北京：科学出版社.

Application of Elastic Foundation Timoshenko Beam Element in ABAQUS

YAN G Zhao1，XU Jian-cong1，YU Jun2
（1.Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education，Tongji University，Shanghai 200092，China；2.School of Civil Engineering and Architecture，Central South University，Changsha 410075，China）

Using the user defined element interface in ABAQUS，the elastic foundation Timoshenko beam element is developed with Fortran language.Comparing with classic analytical solution，the result indicates that the accuracy of the element is high enough.That element not only considers the separation of the foundation spring under tension，but also considers the internal node of curve beam element is not on the line of two end nodes.Using the elastic foundation beam element，the Qingcaosha river-cross water diversion project is analyzed，the results agrees with the measured data of similar projects.

elastic foundation；Timoshenko beam；ABAQUS；element subroutine

TU 471.2

A

1000-5013（2010）04-0448-05

（责任编辑：黄仲一 英文审校：方德平）

2009-09-19

杨钊（1984-），男，博士研究生，主要从事盾构隧道数值计算和模型实验的研究.E-mail：yangzhaolp@126.com.

国家自然科学基金资助项目（40872179）；中国博士后科研基金资助项目（20080440652）