如何在OpenSees中建立墩柱纤维模型并进行PushOver分析

导读：纤维模型在用于结构弹塑性分析时，能以较低的计算成本获得较高的求解精度，且各纤维可以采用材料单轴本构关系，从而避免了确定多维本构关系的困难。墩柱拟静力加载试验是结构抗震经常会做的试验，本文以太平洋地震研究中心（PEER）上边一方形柱试件拟静力加载试验为例（如图1），介绍如何在OpenSees中建立墩柱纤维模型并进行PushOver分析，并附上详细命令流。推荐大家关注我上新精品课《OpenSees从入门到精通21讲：获得OpenSees有限元建模和抗震分析能力》。

图1试件尺寸及配筋

一、节点、约束及坐标转换

# 建模基本命令model BasicBuilder -ndm 3 -ndf 6; # 以平行纸面水平向右为x轴，以竖直向上为y轴，以垂直纸面向外为z轴 file mkdir output; # 建立输出文件夹source LibUnits.tcl; # 定义基本单位及各单位间的换算# 定义节点与约束：node 1 0 0 0; # 以桥墩轴线方向为y轴node 2 0 0 0; # 零长度单元node 3 0 $LCol1 0;node 4 0 $LCol 0;fix 1 1 1 1 1 1 1; # 墩底约束equalDOF 1 2 1 3 5; # 约束横向和沿y轴方向的旋转# 坐标转换 set ColTransfTag 1; # 坐标转换编号set ColTransfType PDelta; # 坐标转换类型，PDelta与Linear的区别就在于考虑了二阶P-Delta效应geomTransf $ColTransfType $ColTransfTag 0 0 -1; # 向量（0，0，-1）是全局坐标的（而此时整体坐标系就是x右y上z垂直纸面向外)，其与x'叉乘（右手规则）得y',进而得z'。所以y'是水平向右的，z'是垂直纸面向内的。

图2 单元划分

二、材料本构

采用可以模拟纵筋屈曲和低周疲劳的Reinforcing Steel模拟钢筋材料，采用可以计及混凝土抗拉强度的Concrete02模拟混凝土材料，采用Bond_Slip01模拟钢筋与混凝土之间的粘结滑移。

# 材料编号 set IDconc 1;         # 核心混凝土set IDconcU 2;        # 保护层混凝土set IDreinf 3;        # 考虑屈曲的纵筋set IDreinf1 4;       # 不考虑屈曲的纵筋set IDBond 5;         # 粘结-滑移材料

1、Reinforcing Steel材料模型

Reinforcing Steel^[1]材料模型可以考虑纵筋屈曲和循环加载导致的强度退化效应,其本构关系由四部分组成，如图3-a)。由于Reinforcing Steel模型考虑了钢筋的初始屈服流幅、包辛格效应、应变软化和刚度退化等力学特征，因此能较好的模拟钢筋在循环加载下的滞回曲线，如图3-b)。图3 Reinforcing Steel 材料模型

#先后定义考虑屈曲和不考虑屈曲两种钢筋材料uniaxialMaterial ReinforcingSteel $IDreinf $Fy $Fu $Es $Esh $esh $eult -GABuck   $lsr   $beta   $r   $gama -CMFatigue $Cf $alpha $Cd;uniaxialMaterial ReinforcingSteel $IDreinf1 $Fy $Fu $Es $Esh $esh $eult -GABuck   $lsr_1   $beta   $r_1   $gama_1 -CMFatigue $Cf_1 $alpha $Cd_1;

2、Concrete02材料模型

Concrete02材料模型基于修正的Kent-Park模型^[2]，可以考虑混凝土的抗拉性能和卸载时的剩余刚度，其本构关系如图4所示。箍筋对混凝土的约束作用通过J. B. Mander约束混凝土模型^[3]加以考虑。

图4 Concrete02本构曲线

# 先后定义了核心混凝土和保护层混凝土材料 uniaxialMaterial Concrete02  $IDconc  $fpc  $epsc0   $fpcu  $epsu  $lambda  $ft  $Et ;uniaxialMaterial Concrete02  $IDconcU  $fpc1  $epsc01  $fpcu1  $epsu1  $lambda  $ft1  $Et1 ;

3、Bond-SP01材料模型

当桥墩试件的钢筋在平台底部锚固长度不长时，应变渗透效应会在墩底一定长度范围内导致桥墩纵筋发生粘结滑移变形，故在该区段附加一个单独的零长度单元，其纵筋采用Bond-SP01^[4]材料模型模拟。Bond-SP01材料模型的应力-滑移骨架曲线与滞回曲线如图5所示，其中和分别为钢筋屈服强度和极限强度，为钢筋应力达到屈服强度时钢筋与混凝土界面处的钢筋滑移值，为钢筋应力达到极限强度时钢筋与混凝土界面处的钢筋滑移值。

图5 Bond_SP01材料模型

uniaxialMaterial Bond_SP01 $IDBond $Fy $Sy $Fu $Su $b $R;

三、纤维截面划分

桥墩截面划分为保护层混凝土、核心混凝土和纵向钢筋三种纤维，如下图6：

图6 纤维截面划分

# 截面编号set ColSecTagFiber 1;     # 钢混段截面编号set ColSecTagFiberBond  2;   # 粘结滑移段截面编号注意：最新的OpenSEES要求定义纤维截面时即指定扭转刚度，否则会报错！section Fiber $ColSecTagFiber -GJ [expr 1.7E08*$Pa*$m4] {# 核心混凝土patch rect       $IDconc      $nfY          $nfZ       -$coreY  -$coreZ  $coreY  $coreZ   # 保护层混凝土（划分纤维的时候边界点顺序先后无妨，y、z坐标都是从小到大，从大到小均可，但定义单元的时候得注意顺序）patch rect  $IDconcU  [expr $nfY+4]  2  -$coverY  -$coverZ  $coverY  -$coreZ   patch rect  $IDconcU  [expr $nfY+4]  2  -$coverY  $coreZ  $coverY  $coverZpatch rect  $IDconcU  2  $nfZ  -$coverY  -$coreZ  -$coreY  $coreZpatch rect  $IDconcU  2  $nfZ  $coreY   -$coreZ  $coverY   $coreZ#定义了墩柱绕z轴弯曲转动时，受压排考虑屈曲和绕y轴弯曲不考虑两种钢筋layer  straight  $IDreinf $numbarColZ  $barAreaCol -$coreYBar $coreZBar -$coreYBar -$coreZBarlayer  straight  $IDreinf $numbarColZ  $barAreaCol  $coreYBar $coreZBar  $coreYBar -$coreZBarlayer  straight  $IDreinf1$numbarColY  $barAreaCol -$coreYBar1   -$coreZBar  $coreYBar1-$coreZBarlayer  straight  $IDreinf1$numbarColY  $barAreaCol -$coreYBar1$coreZBar  $coreYBar1 $coreZBar}

粘结滑移段截面纤维划分与钢混段相同，只是采用的钢筋本构不同，不再赘述。

四、单元定义

1、柔度法单元

与刚度法单元相比，柔度法单元在任意状态下单元力平衡方程都不会产生有限元离散误差，且在相同精度要求下，柔度法需要更少的单元和结构总自由度。采用柔度法单元时，整个桥墩可只采用一个单元，单元内使用4个截面积分点；但为防止应变集中，本文采取Calabrese[5]和Kashani[6]的建议增加一个长度为3倍钢筋屈曲长度的纤维梁柱单元，且单元内使用3个截面积分点，如图3。

# 定义非线性柔度法梁柱单元 ------------------------------------element nonlinearBeamColumn    2      2      3     $numIntgrPts1     $ColSecTagFiber     $ColTransfTag   -iter $maxIters $Tol;   element nonlinearBeamColumn    3      3      4     $numIntgrPts2     $ColSecTagFiber     $ColTransfTag   -iter $maxIters $Tol;

2、零长度单元

顾名思义，零长度单元与普通单元的区别就在于两节点的坐标相同，故长度为零。值得指出的是，两个节点的坐标只是初始相同，受力分析过程中也会有相对位移。

#（由全局坐标系的向量（$x1 $x2 $x3）定义局部坐标系的x轴），同理，由全局坐标系的向量（$yp1 $yp2 $yp3）定义局部坐标系y轴 ,element zeroLengthSection  1  1  2  $ColSecTagFiberBond  -orient  0 1 0 1 0 0;  # 通过此处定义后，零长度单元的局部坐标系与柱单元完全相同了

五、PushOver分析

# 施加轴向力 set IDctrlNode 4;                      # 在结点4上加载；set PCol [expr 968*$kN];              #此处按轴压比=0.1施加轴向力;pattern Plain 1 Linear {load $IDctrlNode 0 -$PCol 0 0 0 0 ;};constraints Transformation ;           # 多点约束，不能用Plainnumberer RCM;system BandGeneral;  set Tol 1.0e-6;                       # 精度（绝对误差）set maxNumIter 200;                  # 最大迭代数set printFlag 1;                      # 计算过程在DOS屏幕上输出信息test NormDispIncr  $Tol  maxNumIter  printFlag; algorithm Newton;set NstepGravity 10;  # number of load stepsset DGravity [expr 1./$NstepGravity];integrator LoadControl $DGravity;analysis Static;analyze $NstepGravity;loadConst -time 0.0;  # 保持（重力）荷载不变，拟时间重置为零# 拟静力加载 （篇幅所限，此处仅为示意，实际需循环往复作用多次）pattern Plain 2 Linear {load $IDctrlNode 1 0 0 0 0 0;};integrator DisplacementControl  $IDctrlNode  1  0.001analyze 2000;

六、记录输出文件

采用Recorder命令提取墩顶位移和墩底反力：

recorder Node  -file output/disp4.out  -node  4  -dof  1 2 3  disp;  # 记录4号节点X、Y、Z方向的位移recorder Node  -file output/reaction2.out  -node  2  -dof  1 2 3  reaction;  # 记录2号节点X、Y、Z方向的反力

值得指出的是，需将recorder命令放在分析命令之前，才能记录到分析结果。最终得到的数值模拟与试验结果的滞回曲线比较如图7所示，可见数值模拟结果与试验符合较好。

图7 数值模拟与试验滞回曲线比较

七、OpenSees从入门到精通

本文完整命令流及其他案例视频讲解可至仿真秀官网本人主页，我的原创精品课《OpenSees从入门到精通21讲》查看。笔者针对OpenSEES学习资料匮乏、上手难的问题，本课程逐行详细讲解如何采用OpenSEES建模与计算分析。课程内容包括命令解释、案例演示以及使用OpenSees进行抗震分析的一些感悟，本课程将持续更新，附件里有全部建模代码和课件，也欢迎感兴趣的朋友加入仿真秀OpenSees软件用户交流讨论（在本公众号对话框回复OpenSees扫码进群交流）。

审核编辑：汤梓红