单弹簧联结单元法的改进算法

颜天佑 李同春 赵兰浩

（1.长江勘测规划设计研究有限责任公司，武汉 430010；2.河海大学 水利水电工程学院，南京 210098）

目前，在分析钢筋与混凝土的相互作用时，所采用的有限元模型主要有三种:分离式、组合式和整体式.采用分离式模型时，常用的联结单元有双弹簧单元、无厚度四边形粘结单元［1］、粘结区单元［2］、界面元［3－4］等.通常联结单元在描述钢筋与混凝土的相互作用时，法向刚度系数kv的取值存在一定随意性，取值过大将会给计算带来一些麻烦，而过小的取值会发生单元相互嵌入问题，致使结构分析模型与实际几何形态有一定的偏差［5］.为解决法向刚度kv的取值问题，文献［6］提出了单弹簧联结单元法.单弹簧联结单元法是在实体单元与梁单元的切向设置单弹簧联结单元，模拟钢筋与混凝土切向的相互作用，其切向刚度同样是由钢筋与混凝土之间的粘结力与滑移量的关系确定；通过建立法向自由度约束方程保证钢筋与混凝土的法向变形协调，避开了人为选择法向刚度系数的困难，而为了方便建立法向自由度之间的关系，实体单元在整体坐标系内求解，梁、杆单元在局部坐标系内求解，能够很方便的考虑钢筋对混凝土的销拴效应.

单弹簧联结单元法中钢筋与混凝土之间的切向力增量主要通过切向刚度与位移增量求得，钢筋和混凝土间的切线刚度随滑移量不断变化，当荷载增量较小时，钢筋与混凝土之间的切向力精确性较高，荷载增量较大时，钢筋与混凝土之间的滑移增量较大，以上一步的切线刚度得到的切向力误差也较大，特别是位移滑移量由本构关系曲线的上升段跃到下降段时，单弹簧联结单元法很难得到一个精确值，一旦迭代试算中的滑移量处于下降段，随着切向力的减小，钢筋混凝土之间滑移量总是不断增加，导致计算不收敛.

本文在单弹簧联结单元法的基础上，在钢筋和混凝土之间直接引入粘结滑移本构关系，由钢筋和混凝土之间的全量滑移量直接求解相互作用力，提出了基于钢筋混凝土粘结滑移关系全量曲线的非线性迭代解法，对单弹簧联结单元法进行了改进，改进后方程求解得到的钢筋混凝土滑移量和切向力严格满足粘结滑移关系，使得单弹簧联结单元法基本不受加载方式的影响，提高了单弹簧联结单元法的求解速度和精度.

1 基于混合坐标的单弹簧联结单元法

单弹簧联结单元法为了计算方便，采用混合坐标系统如图1所示，oxyz 为整体坐标系，o′x′y′z′为某一个梁单元局部坐标系，x′布置在梁的轴向.考虑到实际问题中钢筋可能不在一条直线上，在每一个梁单元节点上建立节点局部坐标系o＊x＊y＊z＊，由相邻梁单元的局部坐标系平均得到.各个坐标系之间的坐标转换矩阵定义如下:R为整体坐标系oxyz与梁单元局部坐标系o′x′y′z′的坐标转换矩阵；r为整体坐标系oxyz与每个定点的局部坐标系o＊x＊y＊z＊的坐标转换矩阵.

若记u、K、F，u′、K′、F′，u＊、K＊、F＊分别为整体坐标系xyz、梁单元局部坐标系x′y′z′、节点局部坐

图1 单弹簧单元混合坐标系示意图

标系x＊y＊z＊内的位移向量、刚度矩阵和荷载向量，则3者之间的关系可用坐标转换矩阵建立如下关系:

梁单元局部坐标系下的有限元平衡方程为

利用式（1）、式（2）代入式（3）中可得到节点局部坐标系x＊y＊z＊中的有限元平衡方程:

其中，K＊＝rKrT＝rRTK′RrT，为梁单元在节点局部坐标系x＊y＊z＊中的刚度矩阵.单弹簧联结单元法应用不协调网格协调位移解法［7］的思想，令双节点的法向自由度的位移相等来直接建立钢筋与混凝土在法线方向的联系.

若记f＊为节点局部坐标系x＊y＊z＊内混凝土对钢筋的作用力，由作用力与反作用力的关系可知，在整体坐标系x，y，z内钢筋对混凝土的作用力为rTf＊.可以建立混凝土和钢筋的增量有限元平衡方程:

由于节点的法向自由度已通过强迫相等的方式保证了钢筋与混凝土之间法向变形协调，即混凝土与钢筋之间的法向作用力已变成内力，不需要在平衡方程中出现，则f＊仅仅为混凝土与钢筋之间的切向粘结力.

式中，τ－s为由实验得到的黏结力与滑移量之间的关系曲线；D为钢筋直径；l为钢筋长度；ds为滑移量；

将式（7）代入（5）、（6），并写成矩阵的形式:

上式即为基于混合坐标系的单弹簧联结单元法的有限元方程.

2 基于滑移关系全量曲线的非线性迭代解法

式（8）为增量方程，钢筋和混凝土间的切线刚度和法线刚度随着滑移量在不断的发生变化，所以对这一类问题进行求解时一般采用增量有限元法，钢筋与混凝土之间的切向力增量主要由式（7）通过切向刚度与增量位移求得.

用增量法求解钢筋与混凝土的切向力时，就是用分段的线性解去逼近非线性解.如何将非线性问题分段线性化，针对不同的情况计算方案也有所不同.为了提高求解精度，荷载增量要取得足够小，这又大大增加了计算量.如果计算步长选择过长，求解的精度往往受到影响，无法知道解的准确程度，有时还造成了计算结构不收敛［8］.另外目前采用的钢筋和混凝土之间粘结滑移本构关系曲线中，有一个明显的下降段，此时如采用一般的迭代法均不能收敛.这主要是因为刚度矩阵不是正定的，存在着“负刚度”的问题.针对这一情况，各国学者提出了不少算法，以克服下降段的不稳定现象.如逐步搜索法，虚加刚性弹簧法，位移控制法，强制迭代法，硬化刚度法，但是各有优缺点［1－2，9］.

基于上述分析，本文提出了基于钢筋混凝土粘结滑移关系全量曲线的非线性迭代解法，在钢筋和混凝土之间直接引入粘结滑移本构关系，不仅可以考虑本构关系曲线的下降段，且不受加载方式的影响，使得求解速度和精度大大提高.具体求解过程如下:

假设第n步荷载已经完成，进入第n＋1步荷载步，第i次迭代已完成，进入第i＋1次迭代.将上一次迭代的结果作为当前状态.对各个计算参量进行更新:

钢筋和混凝土的滑移总量

在单弹簧联结法中，粘结滑移本构关系曲线的切向刚度和界面间滑移量增量得到钢筋和混凝土之间的切向作用力

而本文则根据钢筋和混凝土的滑移总量和粘结滑移关系的全量曲线求得精确的界面切向力，则在局部坐标系x＊y＊z＊内混凝土对钢筋的作用力

在方程（11）中，钢筋和混凝土之间的作用力可以显性的求出，用全量曲线得到的精确的切向力对由增量求得的界面切向力进行替换（如图2所示），即式（11）f＊n＋1替换式（10），然后再对方程（8）的右端项进行调整，将增量荷载以钢筋和混凝土之间的全量作用力表示出来，方程（8）变为可消除由于界面间切向力的改变而产生的不平衡力的非线性方程:

图2 钢筋和混凝土之间的切向力求解示意图

方程求解后，可以求得新的位移增量，进而求得新的总体位移场以及钢筋与混凝土之间新的滑移量，从而可以进入下一次迭代过程，一直到总体位移场和钢筋和混凝土之间的粘结力均同时收敛再转向下一荷载步的计算.如果在计算过程每一次迭代过程采用变刚度迭代，收敛速度更快.

方程（12）中，切向刚度只是求解方程的预判作用，不管切向刚度的正负性，右端项中钢筋与混凝土的作用力与滑移量都是严格的对应关系，通过不断迭代，最终能实现钢筋混凝土的滑移量和切向力满足粘结滑移关系.

3 数值算例

与文献［10］相同，本文以拔出试验作为研究对象，并与双弹簧模型以及文献中的组合式模型进行比较，以验证本文方法的计算精度.拉拔试验模型为C20混凝土块，尺寸为400mm×100mm×100mm，试件中心放置1Φ12钢筋，荷载分为五级施加，等级分别为2kN、4kN、6kN、8kN、10kN.计算模型和网格如图3、图4所示.钢筋和混凝土之间的粘结滑移关系采用Hodue提出的公式:

式中，τ为粘结应力（N／mm2）；s为相对滑移量（mm）；fc为混凝土抗压强度（N／mm2）.

经过计算分析后，得到钢筋应力沿锚固深度分布如图5所示，从图中可以看出采用三种模式的计算结果非常接近，由此说明了网格相近时，本文的单弹簧联结单元法的改进算法，与双弹簧模型以及组合式模型具有相同的计算精度.

图5 钢筋应力分布图

为考察基于滑移关系全量曲线的非线性迭代解法的收敛性，对单弹簧联结单元法改进前后的两种方法进行比较:第一种直接根据增量方程（8）求解，第二种根据改进的方程（12）基于滑移关系全量曲线的非线性迭代解法求解.对上述模型中的钢筋施加12kN的荷载，其中施加荷载级数和刚度变化按照不同的方式进行组合，以分析加载过程和刚度变化对两种不同方法收敛性的影响.分析主要采取了4种求解方式:1级加载常刚度迭代，1级加载变刚度迭代，分10级加载变刚度迭代，分100级加载变刚度迭代.最终以结点不平衡力与外荷载的比值小于10－5作为最终收敛的标准.

图6给出了两种方法按照不同方式进行计算得到的钢筋的应力分布情况，图7给出了钢筋与混凝土之间的切向刚度系数分布情况.从图中可以看出:

（1）当对钢筋施加12kN的荷载时，在钢筋与界面的端部出现了负刚度现象，此时钢筋应力下降的梯度较小.

（2）直接采用增量方程（8）求解时，计算结果受加载步长的影响比较大.荷载1级施加时，计算结果不收敛；荷载分10级施加时，计算结果的收敛性和精度都不高；荷载分100级施加时才具有较好的精度.

（3）当采用基于滑移关系全量曲线的非线性迭代解法求解时，计算精度和收敛性都很好.荷载一级施加时，采用常刚度迭代71次得到的结果与分10、100级施加变刚度迭代得到的结果精度相当；当一级施加采用变刚度迭代时，只需迭代6次就可以得到相同的精度，而且能够很好的反映粘结滑移本构关系曲线中的下降段，说明本文的改进算法受加载方式的影响较小，且具有较好的准确性和计算的高效性.

4 结 论

本文主要对单弹簧联结单元法进行了改进，提出了基于钢筋混凝土滑移关系全量曲线的非线性迭代解法，在钢筋和混凝土之间直接引入粘结滑移本构关系，由钢筋和混凝土之间的全量滑移量直接求解相互作用力，方程求解得到的钢筋混凝土滑移量和切向力严格满足粘结滑移关系，使得单弹簧联结单元法基本不受加载方式的影响，提高了单弹簧联结单元法的求解速度和精度，也使得单弹簧单元法具有更广泛的应用空间.

［1］沈聚敏，王传志，江见鲸.钢筋混凝土有限元与板壳极限分析［M］.北京:清华大学出版社，1993.

［2］吕西林，金国芳，吴晓涵.钢筋混凝土结构非线性有限元理论与应用［M］.上海:同济大学出版社，1996.

［3］杜培荣，卓家寿，孟闻远.钢筋混凝土界面问题分析模型［J］.河海大学学报，2007，35（3）:312－314.

［4］高向玲，李 杰.钢筋与混凝土粘结本构关系的数值模拟［J］.计算力学学报，2005，22（1）:73－77.

［5］戚乐磊，赵鸿铁，梁若筠，等.粘结滑移问题的界面应力元模型［J］.河海大学学报，2004，32（2）:188－191.

［6］赵兰浩，李同春，牛志伟，等.基于混合坐标系的单弹簧联结单元法［J］.河海大学学报:自然科学版，2008，36（6）:796－800.

［7］李同春，李 淼，温召旺，等.局部非协调网格在高拱坝应力分析中的应用［J］.河海大学学报:自然科学版，2003，31（1）:12－45.

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

［9］江见鲸.钢筋混凝土结构非线性有限元分析［M］.西安:陕西科学出版社，1994:91－117.

［10］林新志.考虑粘结滑移的组合式单元模型研究与应用［D］.南京:河海大学，2005.