桥梁支座设计中单元计算模型的开发与实例验证

包善发,陈老伍,蔡 岳

(1.中国市政工程西南设计研究总院有限公司，四川 成都 610000；2.中国公路工程咨询集团有限公司，北京 100089；3.长沙理工大学 土木工程学院，湖南 长沙 410004)

我国已进入第5个地震活跃期[1]，各地地震灾害频发，而桥梁结构隔震装置通过增加系统水平柔度分离上部结构与水平向地震动有效控制结构振动，减少乃至避免结构损伤，具有良好的防震减灾功能。目前，应用较为广泛的隔震装置有天然橡胶支座(LNR)、铅芯橡胶支座(LRB)等，但其仍存在着自身阻尼较小，耗能能力不足，塑性变形较大，自复位能力较弱等缺点。随着近年来形状记忆合金(SMA)的发展，与普通橡胶支座相比，应用其独特的超弹性效应及形状记忆效应设计的SMA橡胶支座在智能隔震领域具有突出表现。

国内外学者就SMA隔震橡胶装置，进行了大量的理论分析与试验研究。李忠献[2]设计了一种基于SMA复合橡胶支座的桥梁智能隔震系统，并与普通橡胶支座进行了减隔震效果对比试验，验证了SMA复合橡胶支座隔震的自恢复能力；庄鹏[3]等设计了一种SMA隔震橡胶支座并研究了其力学特性；左晓宝[4]等建立了超弹性SMA复合阻尼器力学计算模型，分析了滑条槽深、摩擦系数、SMA丝预应变等对阻尼器力学特性的影响规律；ALAM[5]对使用了SMA叠层橡胶支座的三跨连续梁进行了非线性动力分析，对比了不同类型支座对主梁及桥墩位移的影响差异；DEZFULI[6]等提出了新型SMA-LRB支座，并给出了该支座竖向力及水平剪力下的本构模型。随着橡胶、SMA材料本构的明确，SMA橡胶支座的有限元分析愈发精准[7]，但现阶段商用软件中尚未提供适合SMA橡胶支座的单元库，精确分析往往需要建立SMA与橡胶材料复杂超弹性本构关系，并建立支座精细化模型进行动力学分析，一方面工程应用领域学习成本较高；另一方面对于整体结构分析往往涉及多个支座，如需一一建立支座精细化模型则增加了模型复杂度，并随之带来巨量时间成本与硬件成本。

为此，本文基于ABAQUS-UEL，综合考虑了SMA叠层橡胶支座水平剪切性能、竖向压缩刚度、非线性拉伸及临界承载力性能，开发了SMA叠层橡胶支座组合单元模型，并对单元有效性与精确性进行了验证，以低成本高效率的方式实现了SMA叠层橡胶支座的精确模拟，可为相关领域有限元分析提供参考。

1 SMA叠层橡胶支座力学性能

1.1 水平力学性能

与普通橡胶支座近似线弹性水平恢复力不同，SMA叠层橡胶支座(图1)水平力——位移滞回效应具有更为强烈的非线性。BOUC[8]等提出了一种能够表述结构滞回现象的光滑模型，广泛应用于橡胶隔震支座模拟，Casciati[9]在Bouc基础上提出了改进正交方向加卸载判别条件的模型，满足了水平方向各向同性的基本假定。故本文采用Casciati模型表述SMA叠层橡胶支座水平双向力学性能，水平方向恢复力[F]与位移[U]满足如下关系：

图1 SMA叠层橡胶支座构造示意图

[F]=Kd,0[U]+Qd[Z]+cb[U]

(1)

式中:[F]=[Fx,Fy]T；[U]=[Ux,Uy]T；滞回位移[Z]=[Zx,Zy]T；Kd,0为SMA叠层橡胶支座初始屈后刚度;Ux、Uy为水平x，y方向的位移;Qd为支座特征强度;cb为粘滞阻尼。

滞回位移[Z]为位移[U]的函数，满足非线性微分方程[9]：

(2)

同时，根据文献[10]，采用双弹簧模型修正竖向压力下支座水平刚度Kd[11]：

(3)

式中:G为橡胶材料的剪切模量;Ab为剪切面积;ta为支座厚度(不计上下联结钢板);P为竖向压力;Pcr0为支座初始临界承载力。

1.2 竖向力学性能

SMA叠层橡胶支座具有LNR、LRB支座的特性，具有竖向拉压刚度异性的特点，通常计算中往往忽略竖向拉伸刚度的影响，为实现精细化模拟，对竖向拉压刚度分别考虑，采用位移Uz的方向对压缩(Uz≤0)、拉伸(Uz>0)情况分别考虑，即：

(4)

1.2.1竖向压缩刚度

采用线弹性模型模拟SMA叠层橡胶支座的竖向压缩刚度，压缩力Fz1与竖向位移Uz满足：

Fz1=KvUz

(5)

SMA叠层橡胶支座初始竖向压缩刚度Kv1,0主要由橡胶层提供，其计算方法与LRB基本一致[12]：

(6)

(7)

式中：Ecb为橡胶的修正压缩弹模;Ec为橡胶的表观弹模;Eb为橡胶体积弹模。

考虑水平位移U对竖向压缩刚度Kv1的影响，按分段线性公式进行修正[13]：

(8)

式中:R为支座截面外径。

1.2.2竖向拉伸刚度

采用Kumar拉伸模型[14]考虑SMA叠层橡胶支座竖向拉伸性能，拉伸力Fz2按式(9)计算。

(9)

式中:Uz为当前拉伸位移;K为拉伸损伤参数，用以表征支座损伤对承载力降低的影响。根据文献[14]，可取K=10，对于特殊支座结构或特殊材料构成的支座，应根据试验结果重新标定K。

2 单元推导与UEL单元开发

第一节中给出了SMA叠层橡胶支座水平、竖直方向刚度及修正算法，并给出了基于Bonc-Wen模型改进的Casciati模型考虑水平方向恢复力，但由于Casciati模型滞回位移与位移的关系以微分形式给出，不能直接迭代计算，且Casciati模型虽然将Bonc-Wen模型推广到三维应力空间，但其在三维空间中是三轴解耦的，为构建考虑水平向恢复力的三维隔震支座单元，本文提出以两节点三自由度单元与四节点平面应力单元组合而成的组合单元对SMA叠层橡胶支座进行模拟。

2.1 两节点三自由度单元

两节点三自由度单元如图2所示，其i、j结点分别具有x、y、z这3个方向位移分量，将每个单元i、j结点共计6个位移分量表示为一个向量，即：

图2 两节点三自由度单元

{δ}e={δiδj}T={uiviwiujvjwj}T

(10)

结点力可表示为：

{F}e={FiFj}T={FxiFyiFziFxjFxjFxj}T

(11)

对于x、y、z这3个方向，由第一节刚度计算值，结合广义胡克定律，有：

(12)

其中，k1、k2、k3分别为x、y、z这3个方向的刚度，写为矩阵形式有：

{F}e=[k]e{δ}e

(13)

(14)

2.2 多轴Bouc-Wen四节点平面应力单元

2.2.1广义多轴Bonc-Wen模型

文献[15]基于Casciati模型，采用平面应力J2理论，推导出适用于平面应力空间随动强化、等向强化以及混合强化准则的广义多轴Bouc-Wen模型本构关系：

(15)

式中:[σ]为应力张量;[D]为弹性矩阵，二者满足[σ]=[D],[εe]，其中[εe]为弹性应变张量，且[εe]=[ε]-[εp]，[εp]为塑性应变张量；[I]为单位阵；H1、H2为弹性过渡塑性光滑函数，表达式为：

(16)

(17)

式中:φ为J2屈服函数；φ0为J2屈服函数最大值。

[R]为映射矩阵，表达式为：

(18)

现利用该本构关系，推导基于广义多轴Bouc-Wen模型的四节点平面应力单元。

2.2.2四节点平面应力单元推导

对于实际工程结构中不规则曲面上的单元，常采用不规则单元逼近真实解，但对于不规则单元计算较为困难，对此，本文给出形状规则矩形母单元推导，对于不规则单元采用坐标变换方法进行计算，将母单元通过雅克比矩阵J映射于子单元上，本文不再赘述。母单元如图3所示，x方向边长为2a，y方向边长为2b，坐标原点位于矩形形心，其1、2、3、4结点分别具有x、y方向位移分量，即：

图3 平面应力四节点单元

{δ}e={δ1δ2δ3δ4}T=

{u1v1u2v2u3v3u4v4}T

(19)

由文献[16]，设x、y为单元位移，则四节点平面应力单元位移{r}e与形函数[N]分别为：

{r}e=[N]{δ}e

(20)

(21)

(22)

单元应变{ε}可表示为：

{ε}=[B]{δ}e

(23)

[B]=[B1B2B3B4]

(24)

(i=1,2,3,4)

(25)

由虚功原理：

(26)

其中:dV为母单元体积；{F}e为结点力向量；节点虚位移{δ*}T为任意的，故有：

(27)

同时，其微分形式也成立：

(28)

由2.2.1可知平面应力空间随动强化、等向强化以及混合强化准则的广义多轴Bouc-Wen模型本构关系，将式(15)代入式(28)，单元弹性刚度矩阵[Ke]与滞回刚度矩阵[Km]：

(29)

整体刚度矩阵为应力的函数，可写为：

[KP]=α[Ke]+(1-α)[Km]=

(30)

2.3 组合单元

将2.1节推导两节点三自由度单元与2.2节基于广义多轴Bouc-Wen模型的四节点平面应力单元组合，组合单元在局部坐标系x、y、z下位移可表示为：

{δa}={δa1δa2δa3δa4}T

(31)

{δai}={uaivaiwai}T,(i=1,2,3,4)

(32)

由于位移模式包含12个自由度，基于已推导的单元刚度矩阵[k]e及[Kp]，分别对其扩充为12×12矩阵，有：

(33)

(34)

式中:Ki,j与式(30)对应下标Ki,j相同。

则组合单元刚度矩阵为：

2.4 UEL单元开发

基于ABAQUS-UEL可自定义上述组合单元，在增量步开始时主程序将单元结点的外荷载与位移传入UEL子程序，UEL子程序根据传入参数通过刚度矩阵AMATRX更新变量进行相应计算，利用SVARS保存临时状态变量，用RHS保存单元内部结点力，利用RHS与外荷载叠加判断程序收敛性，如收敛则将计算结果传回主程序，如不收敛则主程序进行下一次迭代直至收敛[17]。对于上述组合单元，计算结果的值依赖于内部状态变量(表1)与自由度增量。UEL主要分析流程见图4。

图4 UEL主要分析流程

表1 内部状态变量Table1 Internalstatevariables变量名称变量含义SVARS(1)-SVARS(3)结点1三个自由度方向的残余力SVARS(4)-SVARS(6)结点2三个自由度方向的残余力SVARS(7)-SVARS(9)结点3三个自由度方向的残余力SVARS(10)-SVARS(12)三个剪切方向的变形

3 单元验证

为验证上述单元的有效性与准确性，现根据SMA叠层橡胶支座相关研究成果，对单元水平性能、竖向性能加以验证。

3.1 水平力学性能验证

文献[18]对SMA叠层橡胶支座水平向恢复力进行了试验及理论计算，其中支座主要性能及材料参数:SMA相变温度Mf为-64 ℃，Mx为-40 ℃，Ax为-29 ℃，Af为-5 ℃，SMA弹性模量为4.27×104 MPa;SMA屈服应力为480 MPa;橡胶外径为150 mm;橡胶层厚为5 mm;橡胶层数为7;保护橡胶厚度为5 mm;中间钢板厚度为2.5 mm，橡胶邵氏硬度为45。基于第1节水平恢复力模型，取δ=0.1,β=0.9，对如下两类工况加以验证：

a.工况一：竖向压力荷载P=40 kN，加载频率f=0.02 Hz，支座特征参数屈服力为2.972 9 kN，屈服位移为2.693 3 mm，刚度比为0.488 2。

b.工况二：竖向压力荷载P=60 kN，加载频率f=0.5 Hz，支座特征参数屈服力为3.790 1 kN，屈服位移为2.385 0 mm，刚度比为0.425 3。

文献[18]试验值、模拟值与本文开发单元模拟值在两类工况下的滞回曲线对比见图5。由滞回曲线可计算SMA叠层橡胶支座的等效刚度、单位循环耗能及等效阻尼比，文献[18]试验值、模拟值与本文模拟值见表3。

由图5可知，本文开发单元能良好模拟SMA叠层橡胶支座恢复力-位移的滞回特性，与文献模拟值基本一致，与试验值较为接近，体现了本文开发单元对水平向恢复力的模拟的精确性；由表3可知，分别取不同位移幅值D=5、7、10 mm，本文单元模拟SMA叠层橡胶支座的等效刚度、单位循环耗能、等效阻尼比与文献模拟值、试验值基本一致，等效刚度、单位循环耗能、等效阻尼比与试验值最大误差仅为7%、10%、6%，体现了本文单元对支座水平力学性能模拟的精确性。

(a)工况一对比

表3 等效刚度、单位循环耗能、等效阻尼比对比Table3 Comparisonofequivalentstiffness,unitcycleenergyconsumptionandequivalentdampingratio不同位移幅度值(mm)的等效刚度/(kN·mm-1)D=5D=7D=10文献试验文献模拟本文模拟本文/试验文献试验文献模拟本文模拟本文/试验文献试验文献模拟本文模拟本文/试验工况一0.880.910.901.030.790.860.841.070.710.720.721.02工况二0.971.031.001.030.870.930.911.040.810.770.811.00不同位移幅度值(mm)的单位循环耗能/(kN·mm)D=5D=7D=10文献试验文献模拟本文模拟本文/试验文献试验文献模拟本文模拟本文/试验文献试验文献模拟本文模拟本文/试验工况一13.2613.5913.991.0623.2826.6024.551.0539.0443.1143.011.10工况二17.4719.0418.881.0829.6032.4831.411.0650.1954.6253.641.07不同位移幅度值(mm)的等效阻尼比/%D=5D=7D=10文献试验文献模拟本文模拟本文/试验文献试验文献模拟本文模拟本文/试验文献试验文献模拟本文模拟本文/试验工况一9.609.489.230.969.6310.0410.251.069.919.479.871.00工况二11.4511.7411.661.0211.0411.3911.231.0211.0211.3411.091.01

3.2 竖向力学性能验证

SMA叠层橡胶支座与LNR、LRB支座的主要构造区别在于SMA拉索，而SMA拉索对支座竖向刚度的影响较小，本文对SMA叠层橡胶支座单元竖向刚度上的计算方法与LNR支座一致，众多相关研究如文献[10,12,19]已验证第1节中支座竖向压缩刚度计算方法的正确性，故此处不再重复验证本文开发单元在竖向压缩刚度上的准确性，仅对其拉伸性能进行验证。

根据文献[20]对LRB及RB橡胶支座进行了循环拉伸试验，其支座主要性能及材料参数：橡胶外径为300 mm;橡胶中孔直径为60 mm;橡胶层总厚为57 mm;第一形状系数为25;第二形状系数为4.21;支座总高为120 mm；橡胶剪切模量为0.5 MPa，采用不规则加载方式，加载方案见图6，试验与本文开发单元计算结果对比见图7。

图6 文献[20]拉伸循环加载方案

图7 不规则加载力-位移曲线对比

由图7可知，支座经历反复拉伸时，拉伸位移随着拉伸次数增大而增大，支座拉伸刚度不断降低，这与试验值保持一致，但在卸载段与压缩段衔接部分与试验值存在一定偏差，考虑到试验条件的限制，上述偏差可能是由于无法有效控制加卸载频率导致，总体而言，本文开发单元对于支座竖向拉伸刚度的模拟与试验结果保持一致，具有一定精度。

综上，本文开发组合单元能够良好模拟SMA叠层橡胶支座水平向恢复力以及各向刚度，具有良好性能，相应计算结果与试验值较为接近，在实际有限元计算中可利用本文开发单元对相应支座进行简化，减少不必要的工作量。

4 结语

采用有限元直接建模的方式对SMA叠层橡胶支座进行数值计算较为复杂，而通用有限元软件单元库往往不包含适用于SMA叠层橡胶支座的单元，为满足实际工作的需要，本文开发了SMA叠层橡胶支座单元并进行了验证，得出以下结论：

a.基于SMA叠层橡胶支座水平恢复力以及各向刚度计算方法，推导了两节点三自由度单元与四节点平面应力单元以及二者组合单元刚度矩阵。

b.基于ABAQUS-UEL编制了组合单元子程序，并对该组合单元水平、竖向力学性能计算结果进行了验证，给出了适用于SMA叠层橡胶支座组合单元。对于水平向力学性能，两类工况下滞回曲线模拟值与试验值基本吻合，等效刚度、单位循环耗能、等效阻尼比与试验值最大误差仅为7%、10%、6%；对于竖向力学性能，其反复拉伸力——位移曲线也与试验值基本一致。相应算例表明，使用该单元对SMA叠层橡胶支座进行简化，可以减少建模的难度，降低相应计算成本，且计算结果与相关试验基本一致，可为相关结构计算提供帮助。