基于剪切效应纤维梁单元的结构非线性有限元数值模拟*

李嘉钰，陈梦成，王开心

（华东交通大学 土木建筑学院,南昌 330013）

引 言

基于梁单元的有限元分析模型是近年来人们关注的热点，使用不同梁单元理论的有限元模型对所分析结构的计算精度和计算效率有显著影响.早期的Euler-Bernoulli 经典梁理论[1]应用广泛，在一般情况下可以得到比较满意的结果，但由于未考虑剪切变形给截面带来的影响，在梁高度和跨度相差不大的情况下，会有较大的误差.为了考虑剪切带来的影响，Timoshenko[2]于1921年提出了考虑剪切效应的修正梁理论，该理论仍保留了平截面假定，但实际上梁截面的剪应力与剪应变呈抛物线分布，是不均匀的，这与平截面假定相矛盾，于是引入了截面修正系数，将截面上不均匀的剪应力与剪应变等效为均匀分布，这样既能较好地考虑剪切的影响同时又简化了计算模型.Mari 和Scordelis[3]提出了基于有限元刚度法的纤维梁模型，该模型基于Euler-Bernoulli 梁理论，以单元积分点截面处的力学性能代表整个单元的力学性能，使用线性插值构造积分点的轴向位移场，用Hermite 插值函数构造横向位移场.胡郑州等[4]在UL 列式下，根据连续性介质力学和虚位移原理，并结合Timoshenko 梁理论，推导出了考虑剪切效应的三维纤维梁单元模型.顾观东[5]在Euler-Bernoulli 梁单元模型的基础上，修正了塑性阶段的平截面假定，并引入剪切位移，得到了更精确的梁单元模型.郑欣怡[6]将Timoshenko 梁理论引入单轴本构的弯剪耦合纤维模型中，这种改进后的纤维梁单元模型考虑了弯剪耦合效应和剪切效应对轴向变形的间接影响.Thai 等[7]利用纤维梁单元，对单向轴力和弯矩作用下的桥墩进行了双非线性有限元分析.张传超等[8]基于OPENSEES 软件，对纤维梁单元在钢管混凝土柱应用方面的建模方法进行了研究.马银[9]基于纤维梁单元，对型钢混凝土结构进行建模分析，并进一步推导了考虑黏结滑移关系的型钢混凝土梁柱单元模型.蔺鹏臻等[10]基于ABAQUS 软件，采用精细化的纤维梁单元模型，开发了钢筋和混凝土纤维梁单元材料用户子程序.由此可见，纤维梁单元在理论研究方面和工程应用方面都得到了广泛的应用[11-13].

目前，基于Euler-Bernoulli 梁理论的纤维梁单元其纤维只能考虑单轴受力状态，无法考虑剪切效应；基于考虑剪切效应的Timoshenko 梁理论，多以使用材料多维本构模型来考虑材料非线性，需计算的参数较多，且运算复杂.为此，本文基于考虑剪切效应的三维纤维梁单元有限元模型，推导了该纤维梁单元的截面刚度矩阵、单元弹性刚度矩阵和单元几何刚度矩阵，同时补充了文献[4]在推导插值函数矩阵时，矩阵元素N2,6的遗漏，并根据UL 列式法和弹塑性增量理论，考虑结构的几何非线性和材料非线性，旨在建立压弯剪复杂应力状态下结构非线性有限元分析理论，能在合理描述压弯剪作用的同时得到准确的计算结果.最后，使用MATLAB 编制了相关计算程序，结合钢筋混凝土和矩形钢管混凝土的典型压弯剪构件进行有限元数值模拟，验证了本文所建立理论的有效性、正确性和通用性.

1 考虑剪切效应的纤维模型

1.1 基本假设

本文基于经典纤维模型，并结合Timoshenko 梁理论，建立了考虑剪切效应的纤维模型.纤维梁单元模型如图1所示，并做如下假设：

图1 纤维梁单元Fig.1 The fiber beam element

1） 满足平截面假定，但截面不再与梁中轴线垂直，考虑剪切效应导致的截面翘曲；

2） 钢筋与混凝土充分黏结，无相对滑移，变形协调；

3） 剪应力与剪应变均匀分布.

1.2 刚度矩阵

1.2.1 截面刚度矩阵

对于截面内任一纤维单元，其空间几何变换矩阵为

与其材料性质相关的本构关系矩阵为

其中，Ei和Gi分别为第i个纤维单元的弹性模量和剪切模量.

通过积分运算，得到任一纤维单元的刚度：

其中，Ai为纤维单元的截面面积.

对截面内的全部纤维单元进行积分，得到截面刚度矩阵：

将式(1)、(2)代入式(4)中，展开并进行积分运算，可以得到截面刚度矩阵的显式为

1.2.2 弹性刚度矩阵

在得到截面刚度矩阵的基础上，可采用位移形函数将结构的位移与截面的变形联系起来，从而推导三维纤维梁单元的刚度矩阵.每个梁单元共有2 个节点，每个节点有6 个自由度，单元节点位移向量为

其中，(u1,u2),(v1,v2),(w1,w2)分别为节点在x轴、y轴和z轴方向的线位移，(θ1x,θ2x)为节点的扭转角位移，(θ1y,θ2y)和(θ1z,θ2z)分别为节点在xOz和xOy平面内的转动角位移.

单元轴线上任一点x处截面的位移用单元两端节点位移可表示为

式中，N(x)为单元形函数矩阵.轴向位移及扭转角位移均采用相同的线性位移模式，横向弯曲位移均采用三次式的位移模式[14]，此时单元位移插值函数可写为

单元形函数矩阵为

其中

对单元形函数矩阵进一步推导，可得到线性应变矩阵BL，再由虚功原理和变分原理得单元弹性刚度矩阵为

做积分运算后，得到单元弹性刚度矩阵的显式为

其中，单元弹性刚度矩阵为对称矩阵，即

1.2.3 几何刚度矩阵

在两节点12 个自由度单元线弹性分析力学模型的基础上，本文将建立单元的几何非线性分析力学模型，并推导出单元的几何刚度矩阵.

非线性有限元方程为

其中，BNL为非线性应变矩阵，ΔP为等效荷载.

在式(12)中，左边第一项积分为弹性刚度矩阵，通过第二项积分，可推导几何刚度矩阵.

其中，ATσ=Mθ=MGΔue，A为非线性应变与位移梯度矢量之间的转换矩阵，M为应力分量矩阵，θ为位移梯度矢量，为位移梯度矢量与单元节点位移矢量之间的转换矩阵，故有G

因此，单元几何刚度矩阵为

式中G和的显式为[4]

其中，σ11为沿x轴的轴向应力，τxy为垂直x轴沿y轴的剪应力，τxz为垂直x轴沿z轴的剪应力，τyz为垂直y轴沿z轴的剪应力.

通过式(15)的积分运算可得单元几何刚度矩阵，最终，单元刚度矩阵为弹性刚度矩阵和几何刚度矩阵之和，即

2 复杂应力状态下纤维单元应力的更新

2.1 等效应力与等效应变

本文将采用Mises 屈服准则判断复杂应力状态下纤维单元的屈服状态，并据此理论求等效应力与等效应变.式(17)和(18)为等效应力与等效应变的分量表示方法：

等效应力与等效应变存在以下关系：

其中，φ为已知的材料本构关系.

2.2 基于增量理论应力的更新

由于塑性本构关系中的应力和应变不存在一一对应的关系，一般只能建立应力与应变增量之间的关系，以增量形式来表示.依据此理论对纤维单元应力状态的更新步骤如下：

1） 根据已知的应力分量 σ和应变分量ε 计算出等效应力和 等效应变，计算应变分量增量 dε.

2） 在Mises 屈服准则下，通过比较等效应力与纤维单元的屈服应力判断是否屈服，若未屈服，按第3）步弹性理论更新应力；若屈服，按第4）步塑性理论更新应力.

3） 根据弹性理论，计算弹性Jacobi 矩阵De：

其中，E为纤维单元的弹性模量，µ为Poisson 比.

然后，计算应力分量增量 dσ：

最后，按第5）步更新应力.

4） 根据塑性理论，计算塑性Jacobi 矩阵Dep：

其中

式(23)中，G为纤维单元的剪切模量；为通过式(17)所求得的等效应力；Ep为塑性模量，它与弹性模量E和切线模量Et之间的关系为Ep=EEt/(E−Et)，切线模量可通过对本构关系的求导得到.

平均应力、应力偏量按式(24)计算：

然后，计算应力分量增量 dσ：

最后，按第5）步更新应力.

5） 得到纤维单元的应力增量后，更新应力：

更新材料应力状态的流程图如图2所示.

图2 增量理论应力更新流程图Fig.2 The flow chart for stress updating with the incremental theory

3 算例分析

3.1 算例一

一正方形截面钢筋混凝土柱[15]，如图3所示，柱高1.65 m，截面尺寸为550 mm×550 mm，纵向受力钢筋直径为20 mm，等距分布，保护层厚度为40 mm，考虑结构自重的影响，将重力等效为柱顶的集中力P，同时，在顶部作用一水平集中荷载V，考虑几何非线性和材料非线性，对此结构做水平力-侧移的非线性全过程分析.

图3 钢筋混凝土柱（单位：mm）Fig.3 The reinforced concrete column (unit:mm)

算例分析采用本文所给出的三维纤维梁单元建模计算，计算时将原结构划分为10 个单元，单元的每个节点有6 个自由度；划分纤维单元时，尽可能保证原钢筋面积与纤维网格的面积相等，优先考虑此条件；将截面划分为31×31的正方形纤维网格，网格按照结构的实际情况分为钢筋纤维和混凝土纤维，如图4所示，其余参数见表1.

表1 计算参数Table 1 Calculation parameters

图4 纤维单元截面Fig.4 The fiber element section

本文迭代求解的方法采用杜修力等提出的位移控制新方法[16]，控制点为水平力施加点，位移增量控制步步长=1 mm，总位移增量步数n=25，迭代收敛准则采用位移收敛准则，允许误差δ=10−3.

3.1.1 混凝土本构模型

考虑到混凝土受到箍筋约束作用，故在数值模拟中，受拉区混凝土采用经Scott 等[17]修正后的Kent-Park[18]单轴混凝土本构模型，如图5所示.需要说明的是，对于受拉区混凝土，经试算表明，可以忽略其受拉贡献，即取受拉混凝土纤维的拉应力为0，这一简化在不影响计算精度的前提下提高了计算效率.混凝土抗压强度及相应的混凝土峰值应变和极限应变均取文献[15]中所给出的数值.混凝土受压本构模型具体表达式如下：

图5 混凝土本构模型Fig.5 The constitutive model of concrete

其中

式中，K为考虑箍筋约束效应所引起的混凝土强度增大系数，根据本算例参数可计算得K=1.21；ρsv为箍筋的体积配筋率；fyh为箍筋的屈服强度；εu为混凝土极限压应变，本文取 εu=0.024 8；Zm为约束混凝土应力应变曲线下降段的斜率，经式(29)计算，本文取Zm=13.66.

3.1.2 钢筋本构模型

钢筋的本构模型采用双折线弹塑性本构模型，如图6所示，此本构模型具体表达式如下：

图6 钢筋双折线本构模型Fig.6 The constitutive model of steel reinforcement

式中，Es为钢筋的弹性模量；fy为钢筋的屈服强度；εy为钢筋的屈服应变；k为钢筋硬化段斜率，本文与文献[15]相同，取k=0.01Es.

3.1.3 数值模拟结果

利用弹塑性增量理论结合考虑剪切效应下的三维纤维梁单元模型，编制相应的MATLAB 程序计算，将数值模拟结果与文献[15]的数据进行对比分析，如图7和表2所示.

图7 荷载-位移关系曲线Fig.7 The load-displacement curves

表2 数值模拟结果Table 2 Numerical simulation results

由图7和表2可见：本文的计算结果与文献[15]结果吻合较好，水平荷载的全过程有限元数值模拟结果略低于文献[15]的结果，其相对误差均在8%以内，结构的弹性阶段、弹塑性阶段以及下降段均较为明显.上述结果说明，本文所建立的复杂应力状态下考虑剪切效应的纤维梁单元理论是正确的，依据该理论所编制的非线性有限元分析程序，能有效地解决实际问题.

3.2 算例二

一矩形钢管混凝土压弯剪构件[19]，试件长度L=1 200 mm，截面高D=300 mm，截面宽B=200 mm，外钢管厚度t=5.6 mm，含钢率 α=0.1，所用混凝土强度等级为C60，钢管为Q345 钢，试件端部所受轴压N=0.4Nu（Nu为该构件的轴压极限承载力），跨中强轴方向受水平荷载V.在数值模拟有限元分析模型中，试件两端均为铰接，钢材的本构模型采用双折线模型，混凝土的本构模型采用塑性损伤模型[20].计算简图如图8所示，采用本文所给出的考虑剪切效应的纤维梁单元模型建模，使用位移控制新方法[16]进行有限元非线性方程的迭代求解，计算水平荷载V与支座处转角位移R(R=2Δ/L，Δ为外荷载V产生的水平位移）的全过程关系曲线.

图8 矩形钢管混凝土压弯剪试件计算模型Fig.8 The calculation model for the rectangular CFST column under the compression-bending-shear loading condition

图9为本文数值模拟结果与文献[19]的V-R曲线.由图可知，两者在弹性阶段吻合程度良好，在弹塑性阶段和下降段，由于纤维梁单元模型不能模拟钢管和混凝土的相互作用，本文的计算结果稍滞后于文献[19]的结果，但偏差相对较小，从总体上来说，两者结果较为吻合.此外，本文计算出的构件所受最大横向力为479.8 kN，文献[19]的最大横向力约为480.1 kN，说明本文的纤维梁单元模型能较好地模拟构件的横向极限承载力.

图9 矩形钢管混凝土试件V-R 曲线Fig.9 V-R curves of the rectangular CFST column

综上所述，本文所建立的复杂应力状态下结构非线性分析理论是正确的，证明了依据此理论所编制的非线性程序是可行、有效的.

4 结 论

1） 算例分析的结果表明，本文基于考虑剪切效应的纤维梁单元，所建立的压弯剪复杂应力状态下结构非线性有限元分析理论是通用、可行和正确的.

2） 利用本文理论进行结构的非线性全过程分析时，可较好地模拟出结构的弹性阶段、弹塑性阶段和下降段，且特征点较为明显.

3） 依据本文理论所编制的非线性有限元分析程序是有效、可靠的，可以较好地解决理论研究和工程应用中所涉及的结构非线性问题.