近年来, 随着我国科技进步和社会发展, 各种高楼和大跨建筑拔地而起, 建筑形式和体量的发展到达了前所未有的高度。然而随之而来的是, 越来越多的结构需要进行罕遇地震作用下的弹塑性分析以确定其安全性。由于强震作用下, 结构材料往往会进入塑性阶段, 若仍然采用弹性假定进行分析, 则计算结果会与结构实际响应差别很大。因此, 需要对结构的非线性性能进行合理分析, 才能较为准确地得到结构的大震性能指标。

　　 纤维梁单元是建筑结构领域常用的非线性单元, 其将梁单元截面分割成一定数量的纤维, 每条纤维仅可进行轴向变形, 并且可以采用不同的单轴材料本构。同时, 截面本身满足平截面假定, 截面刚度通过对截面上每个纤维进行积分得到。这种单元既可以较为准确地模拟梁单元的非线性性能, 同时计算量也适中, 因此广泛应用于建筑结构的梁柱单元非线性模拟中^[1,2,3]。然而, 在实际使用纤维单元时, 有时会出现计算不准确、结果不合理的情况, 这大部分是由于纤维梁单元选用不合理导致的。

　　 本文推导了几种纤维梁单元的计算原理并研究了具体实现方法, 这些单元包括不考虑剪切变形的欧拉梁以及考虑剪切变形的铁木辛柯梁, 其中铁木辛柯梁又可以按照结点数量和插值函数分为2, 3, 4结点单元。文中分别给出了各类单元的数值算例, 并结合ABAQUS软件进行了分析对比, 给出了纤维梁单元选择的一般性建议。目前, 这些单元已在北京市建筑设计研究院有限公司自主研发的设计软件Paco中完整实现, 验证结果令人满意。

　　 假设梁单元的两个结点为1和2, 分别包含6个自由度。单元坐标系原点取为结点1, 单元坐标系x轴指向结点2。则其结点位移向量u_i定义如下:

　　 $\mathit{u}{}_i=\left\{\begin{array}{cccccc}u{}_i& v{}_i& w{}_i& \theta {}_{\text{x}i}& \theta {}_{\text{y}i}& \theta {}_{\text{z}i}\end{array}\right\}{}_{i=1,2}^{\text{Τ}}(1)$

　　 式中:u_i, v_i, w_i分别为i号结点的三个平动位移;θ_xi, θ_yi和θ_zi分别为i号结点的三个转动位移。

　　 则单元位移向量u_e为:

　　 $\mathit{u}{}_{\text{e}}=\left\{\begin{array}{c}\mathit{u}{}_1\\ \mathit{u}{}_2\end{array}\right\}(2)$

　　 梁截面位移向量u (x) 可写为:

　　 $\mathit{u}(x)=\left\{\begin{array}{cccccc}u& v& w& \theta {}_{\text{x}}& \theta {}_{\text{y}}& \theta {}_{\text{z}}\end{array}\right\}{}^{\text{Τ}}(3)$

　　 式中:x为计算截面位置到结点1的距离;u, v, w和θ_x, θ_y, θ_z分别为与结点1距离为x处截面的三个平动位移和三个转动位移。

　　 经典欧拉梁单元满足平截面假定, 且截面变形后仍垂直于截面形心线。对于梁截面的轴向位移和扭转位移, 采用线性插值函数, 其形式如下:

　　 $\begin{array}{l}u=[\begin{array}{cc}L{}_1& L{}_2\end{array}]\left\{\begin{array}{cc}u{}_1& u{}_2\end{array}\right\}{}^{\text{Τ}}(4)\\ \theta {}_{\text{x}}=[\begin{array}{cc}L{}_1& L{}_2\end{array}]\left\{\begin{array}{cc}\theta {}_{\text{x}1}& \theta {}_{\text{x}2}\end{array}\right\}{}^{\text{Τ}}(5)\\ \{\begin{array}{c}L{}_1(x)=x/L\hfill \\ L{}_2(x)=1-x/L\hfill \end{array}(6)\end{array}$

　　 式中:L₁和L₂为线性插值函数;L为梁长。

　　 对于梁截面的挠度采用三次Hermite插值函数, 其形式如下:

　　 $\begin{array}{l}v=[\begin{array}{cccc}{\rm H}{}_1& {\rm H}{}_2& {\rm H}{}_3& {\rm H}{}_4\end{array}]\left\{\begin{array}{cccc}v{}_1& \theta {}_{\text{z}1}& v{}_2& \theta {}_{\text{z}2}\end{array}\right\}{}^{\text{Τ}}(7)\\ w=[\begin{array}{cccc}{\rm H}{}_1& -{\rm H}{}_2& {\rm H}{}_3& -{\rm H}{}_4\end{array}]\left\{\begin{array}{cccc}w{}_1& \theta {}_{\text{y}1}& w{}_2& \theta {}_{\text{y}2}\end{array}\right\}{}^{\text{Τ}}(8)\end{array}$

　　 其系数表达式如下:

　　 $\{\begin{array}{c}{\rm H}{}_1(x)=1-3x{}^2/L{}^2+2x{}^3/L{}^3\hfill \\ {\rm H}{}_2(x)=x-2x{}^2/L+x{}^3/L{}^2\hfill \\ {\rm H}{}_3(x)=3x{}^2/L{}^2-2x{}^3/L{}^3\hfill \\ {\rm H}{}_4(x)=-x{}^2/L+x{}^3/L{}^2\hfill \end{array}(9)$

　　 式中H_i (i=1, 2, 3, 4) 为三次Hermite插值函数。

　　 根据平截面假定, 梁截面任一点轴向变形$\overline{u}$为:

　　 $\overline{u}(x,y,z)=u+z\theta {}_{\text{y}}-y\theta {}_{\text{z}}(10)$

　　 则截面任一点正应变ε_xx为:

　　 $\epsilon {}_{\text{x}\text{x}}=\frac{\partial \overline{u}(x,y,z)}{\partial x}=u^{\prime }+z\theta {}_{{\text{y}}^{\prime }}-y\theta {}_{{\text{z}}^{\prime }}(11)$

　　 假设截面扭转为弹性并均匀分布, 且扭转变形不与其他变形相耦合, 引入截面扭率α=θ_x′, 则应变向量ε可写为:

　　 $\begin{array}{l}\mathit{\epsilon }=\left\{\begin{array}{cc}\epsilon {}_{\text{x}\text{x}}& \alpha \end{array}\right\}{}^{\text{Τ}}=\mathit{B}\mathit{u}{}_{\text{e}}(12)\\ \mathit{B}=\left[\begin{array}{cc}\mathit{B}{}_1& \mathit{B}{}_2\end{array}\right](13)\\ \mathit{B}{}_i=\left[\begin{array}{cccccc}L{}_{i^{\prime }}\hfill & -y{\rm H}{}_{2i-1}˝\hfill & -z{\rm H}{}_{2i-1}˝\hfill & 0\hfill & z{\rm H}{}_{2i}˝\hfill & -y{\rm H}{}_{2i}˝\hfill \\ 0\hfill & 0\hfill & 0\hfill & \frac{J}{A}L{}_{i^{\prime }}\hfill & 0\hfill & 0\hfill \end{array}\right]{}_{(i=1,2)}(14)\end{array}$

　　 式中:B为单元应变矩阵;B₁和B₂分别为对应于结点1和结点2的应变矩阵;J为截面抗扭惯量;A为截面面积;L_i′为线性插值函数L_i关于x的一阶导数;H_i″为三次Hermite插值函数H_i关于x的二阶导数。

　　 根据虚功原理:

　　 ${\displaystyle {\int }_V\delta }\mathit{\epsilon }{}^{\text{Τ}}\mathit{\sigma }\text{d}V=(\delta \mathit{u}{}_{\text{e}}){}^{\text{Τ}}\mathit{{\rm P}}(15)$

　　 式中:σ为应力向量;P为节点荷载向量;V为单元体积。

　　 将式 (12) 代入式 (15) 中, 并考虑δu_e任意性, 可得:

　　 ${\displaystyle {\int }_V\mathit{B}}{}^{\text{Τ}}\mathit{\sigma }\text{d}V=\mathit{{\rm P}}(16)$

　　 将式 (16) 写为增量形式:

　　 ${\displaystyle {\int }_V\mathit{B}}{}^{\text{Τ}}\text{Δ}\mathit{\sigma }\text{d}V=\text{Δ}\mathit{{\rm P}}(17)$

　　 将本构关系Δσ=DΔε和式 (12) 代入式 (17) , 可得:

　　 $\mathit{k}{}_{\text{e}}\text{Δ}\mathit{u}{}_{\text{e}}=\text{Δ}\mathit{{\rm P}}(18)$

　　 其中k_e为单元刚度矩阵, 表达式为:

　　 $\mathit{k}{}_{\text{e}}={\displaystyle {\int }_V\mathit{B}}{}^{\text{Τ}}\mathit{D}\mathit{B}\text{d}V(19)$

　　 式中D为切线弹性矩阵。

　　 假设截面分为n条纤维, 取每个纤维形心点的本构关系作为整个纤维的本构关系, 则截面的刚度矩阵可写为:

　　 $\mathit{k}{}_{\text{e}}\approx {\displaystyle {\int }_0^L\left({\displaystyle \sum _{i=1}^n\mathit{B}}{}_i^{\text{Τ}}\mathit{D}{}_i\mathit{B}{}_i\text{Δ}A{}_i\right)}\text{d}x(20)$

　　 式中:B_i为纤维i形心处的单元应变矩阵;D_i为纤维i的切线弹性矩阵;ΔA_i为纤维i的面积。

　　 取纤维i的轴向切线刚度为E_ti, 截面剪切刚度为G, 则D_i矩阵可写为:

　　 $\mathit{D}{}_i=[\begin{array}{cc}E{}_{\text{t}i}& \\ & G\end{array}](21)$

　　 对杆件沿长度方向采用Gauss-Legendre积分求解, 即可求出单元刚度矩阵k_e, 从而用于弹塑性分析。

　　 铁木辛柯梁允许梁产生剪切变形, 即梁截面变形后并不一定垂直于形心轴线。对于梁变形的6个自由度, 均采用线性插值。则梁截面位移向量u (x) 定义如下:

　　 $\mathit{u}(x)=\left[\begin{array}{cc}L{}_1\mathit{{\rm I}}& L{}_2\mathit{{\rm I}}\end{array}\right]\mathit{u}{}_{\text{e}}(22)$

　　 式中I为6阶单位矩阵。

　　 根据平截面假定, 梁截面任一点的三向变形为:

　　 $\{\begin{array}{c}\overline{u}(x,y,z)=u+z\theta {}_{\text{y}}-y\theta {}_{\text{z}}\hfill \\ \overline{v}(x,y,z)=v-z\theta {}_{\text{x}}\hfill \\ \overline{w}(x,y,z)=w+y\theta {}_{\text{x}}\hfill \end{array}(23)$

　　 式中$\overline{v}$和$\overline{w}$为截面任一点两个方向的挠度变形。

　　 则应变可写为:

　　 $\{\begin{array}{c}\epsilon {}_{\text{x}\text{x}}=\frac{\partial \overline{u}(x,y,z)}{\partial x}=u^{\prime }+z\theta {}_{{\text{y}}^{\prime }}-y\theta {}_{{\text{z}}^{\prime }}\hfill \\ \gamma {}_{\text{x}\text{y}}=\frac{\partial \overline{v}(x,y,z)}{\partial x}+\frac{\partial \overline{u}(x,y,z)}{\partial y}=v^{\prime }-z\theta {}_{{\text{x}}^{\prime }}-\theta {}_{\text{z}}\hfill \\ \gamma {}_{\text{x}\text{z}}=\frac{\partial \overline{w}(x,y,z)}{\partial x}+\frac{\partial \overline{u}(x,y,z)}{\partial z}=w^{\prime }+y\theta {}_{{\text{x}}^{\prime }}+\theta {}_{\text{y}}\hfill \end{array}(24)$

　　 式中γ_xy和γ_xz为截面任一点y向和z向的剪应变。

　　 假设截面扭转为弹性并均匀分布, 且扭转变形不与其他变形耦合。另假设横向剪切弹性, y向和z向剪切面积分别为A_y和A_z, 则应变向量ε可写为:

　　 $\begin{array}{l}\mathit{\epsilon }=\left\{\begin{array}{cccc}\epsilon {}_{\text{x}\text{x}}& \alpha & \gamma {}_{\text{x}\text{y}}& \gamma {}_{\text{x}\text{z}}\end{array}\right\}{}^{\text{Τ}}=\mathit{B}\mathit{u}{}_{\text{e}}(25)\\ \mathit{B}=\left[\begin{array}{cc}\mathit{B}{}_1& \mathit{B}{}_2\end{array}\right](26)\\ \mathit{B}{}_i=\left[\begin{array}{cccccc}L{}_{i^{\prime }}\hfill & 0\hfill & 0\hfill & 0\hfill & zL{}_{i^{\prime }}\hfill & -yL{}_{i^{\prime }}\hfill \\ 0\hfill & 0\hfill & 0\hfill & \frac{J}{A}L{}_{i^{\prime }}\hfill & 0\hfill & 0\hfill \\ 0\hfill & \frac{A{}_{\text{y}}}{A}L{}_{i^{\prime }}\hfill & 0\hfill & 0\hfill & 0\hfill & -\frac{A{}_{\text{y}}}{A}L{}_i\hfill \\ 0\hfill & 0\hfill & \frac{A{}_z}{A}L{}_{i^{\prime }}\hfill & 0\hfill & \frac{A{}_{\text{y}}}{A}L{}_i\hfill & 0\hfill \end{array}\right]{}_{(i=1,2)}(27)\end{array}$

　　 根据第1.1节的推导, 假设截面分为n条纤维, 单元刚度矩阵k_e可写为:

　　 $\mathit{k}{}_{\text{e}}\approx {\displaystyle {\int }_0^L\left({\displaystyle \sum _{i=1}^n\mathit{B}}{}_i^{\text{Τ}}\mathit{D}{}_i\mathit{B}{}_i\text{Δ}A{}_i\right)}\text{d}x(28)$

　　 取纤维i上的切线刚度E_ti, 并假设剪应力在截面上均匀分布, 剪切刚度为G, 则D_i矩阵可写为:

　　 $\mathit{D}{}_i=\left[\begin{array}{cccc}E{}_{\text{t}i}& & & \\ & G& & \\ & & G& \\ & & & G\end{array}\right](29)$

　　 与欧拉梁相同, 对杆件沿长度方向积分, 可以求出单元刚度矩阵k_e, 从而用于弹塑性分析。

　　 第1节对欧拉梁和2结点铁木辛柯梁进行了公式推导, 此外, 根据插值函数和结点数的不同, 铁木辛柯梁还包括3结点和4结点梁。目前, 这些单元已在北京市建筑设计研究院有限公司自主研发的结构设计软件Paco^[4]中实现, 以下是算例证明。

　　 图1 (a) 所示悬臂柱, 高6m, 截面为工字形, 尺寸如图1 (b) 所示, 材料采用Q235钢材, 弹性模量为206GPa, 泊松比为0.3, 钢材本构采用随动强化模型 (图1 (c) ) , 屈服强度取标准值235MPa, 屈服后刚度比为0.01, 阻尼比取0.02, 密度为0。在柱顶施加平移质量1t。

　　 输入地震动为El Centro波, 原始波形如图1 (d) 所示, 持时12s, 输入步长为0.02s, 为了计算悬臂柱塑性行为, 峰值加速度 (PGA) 调整为2 500cm/s²。动力时程分析采用常加速度Newmark-β法。

　　 采用Paco和ABAQUS分别建模计算, 并根据单元类型不同分组, 组别和特征如表1, 2所示。

　　 由表1, 2可知, ABAQUS的三种空间梁单元在Paco中均有对应实现的单元, 可以进行对比, 此外Paco还提供4结点铁木辛柯单元作为选择。两种软件的位移时程对比结果如图2所示。可以看出, 在单元数划分合理的情况下, Paco与ABAQUS软件在该悬臂梁模型下的计算结果基本吻合, 证明了Paco实现的两类非线性梁单元的正确性。

　　 图3为两种纤维梁单元在杆件不同细分情况下的x向位移时程曲线。从图3可以看出, 欧拉梁细分前和细分后有较大差别 (细分为1段和3段时, 顶点最大位移相差22.9%) , 这与弹性理论的结论不符。这是由于三次Hermit插值函数不能很好地描述底部屈服后单元的曲率分布^[5], 因此在构件端部屈服之后, 采用Hermit函数计算会产生误差。同时, 欧拉梁无法考虑剪切变形的影响, 因此不适用模拟短梁短柱, 建议仅将其应用在细长梁柱的弹塑性分析中, 并应进行细分。

　　 对于铁木辛柯梁, 采用一个单元模拟也不是很准确, 但误差随着单元内结点数量的增多而减少, 将2, 3, 4结点铁木辛柯梁单元细分1段与3段相比, 顶点最大位移差分别为26.63%, 11.8%和5.7%。采用多结点单元可以通过较少的划分次数获得合理结果, 但单元相对复杂、计算速度偏慢, 而采用2结点单元则需要多次划分, 才能保证正确性, 但单元简单、计算速度快。此外, 多结点单元形函数阶次较高, 相对于划分后呈折线形的2结点单元, 能够更好地模拟梁单元的变形形态。因此在工程实际中, 应根据计算精度和效率的要求, 选择较少划分的多结点单元或较多划分的2结点单元。

　　 对于3结点和4结点铁木辛柯梁, 由于属于多结点单元, 存在单元内结点。这些结点与其他单元无关, 只存在于当前单元中。为了减少出口自由度的数量, 常采用静力凝聚的方法, 将这些结点的自由度凝聚到端部结点中。为了验证静力凝聚在本单元上的效果, 采用图1所示模型, 分别对3结点和4结点铁木辛柯梁在有无静力凝聚的情况下进行测试, 测试结果如图4, 5所示。从图中可以看出, 在PGA=250cm/s²时, 悬臂柱处在弹性阶段, 此时是否采用静力凝聚对结点位移时程没有影响, 两条曲线完全重合, 说明静力凝聚在弹性状态下效果良好。而对于PGA=2 500cm/s²的情况, 两条曲线不再重合, 且差别较大。

　　 不仅如此, 对采用静力凝聚的情况下, 不同的单元细分情况进行测试, 结果如图6所示。

　　 可以看出, 若采用静力凝聚方法, 两种多结点单元均不能通过细分单元得到收敛结果。由于静力凝聚会将内部结点的刚度分配到端部结点上, 会产生一定的误差, 而这种误差在结构进入塑性阶段后不可估计, 从而可能导致结果不合理。因此, 在采用3结点或4结点铁木辛柯纤维梁进行弹塑性动力时程分析时, 应慎重采用静力凝聚方法, 而采用内部结点出口的方式则结果良好。

　　 采用纤维梁单元能够较好地模拟梁单元的非线性性能, 然而纤维划分不当的情况下, 可能会使结构响应出现一定的差别。例如, 针对下列工字形截面, 常用的划分方式如图7 (a) 所示, 翼缘和腹板均采用一维划分, 在厚度方向不再划分。采用这种方式, 在翼缘和腹板较薄时, 误差并不大, 但当翼缘或腹板较厚时 (图7 (b) ) , 由于计算时忽略了纤维本身的抗弯能力, 因此计算截面抗弯刚度会偏小, 采用这种划分方式就会产生较大误差。

　　 为了减小这种误差, 可采用两种方法:1) 方法1:将板件沿厚度方向进行划分 (图7 (c) ) 。这种方法可以直接解决上述问题, 但会增加纤维数量, 且较难确定合理的划分层数。2) 方法2:在计算纤维刚度时, 采用其他方法考虑纤维的抗弯贡献, 比如修正纤维位置, 使纤维截面对两个主轴的惯性矩与原截面相同等。这种方法能在不增加单元数量的情况下考虑板厚效应, 但需选择合理的调整方法。

　　 为了验证纤维调整对于结果的影响, 仍然采用图1所示算例, 更改截面尺寸和顶点质量大小建立测试组 (表3) , 地震动PGA=4 000cm/s², 进行弹塑性时程分析, 纤维调整采用方法2分析结果见表4。

　　 从表4可以看出, 随着翼缘腹板变厚, 调整前和调整后的结构响应差别会越来越大, 因此在采用纤维梁单元进行非线性时程分析时, 应注意考虑板厚造成的影响。

　　 对于建筑结构大震弹塑性分析来说, 结构的变形很大, 几何非线性效应明显。因此, 必要时应在采用纤维梁单元时考虑几何非线性。要实现纤维梁单元的几何非线性, 主要有如下几个要点:

　　 (1) 建议采用更新拉格朗日法列式 (U.L.法) , 一般情况下, 其应变矩阵的表达式更为简单, 但在每一步积分过程中需要不断更新其参考构型, 也就是平衡方程的建立位置需要逐步更新。

　　 (2) 平衡方程中需引入非线性项的影响, 其中最主要的非线性项为几何刚度k_NL, 其一般表达式如下:

　　 $\mathit{k}{}_{\text{Ν}\text{L}}={\displaystyle {\int }_V\mathit{B}}{}_{\text{Ν}\text{L}}^{\text{Τ}}\mathit{\sigma }\mathit{B}{}_{\text{Ν}\text{L}}\text{d}V(30)$

　　 式中:B_NL为非线性应变与位移的转换矩阵;σ为Cauchy应力^[6]。

　　 本文推导了常见纤维梁单元的计算方法和公式, 包括不考虑剪切变形的欧拉梁和考虑剪切变形的铁木辛柯梁, 并将这些单元进行了软件实现, 得到了以下结论:

　　 (1) 在进行弹塑性时程分析时, 欧拉纤维梁也需要划分多个单元才能获得准确结果, 且不能考虑剪切变形, 因此仅建议其在模拟细长梁时使用, 且须进行细分;2, 3, 4结点的铁木辛柯梁虽同样需要多段划分, 但多结点单元所需划分次数较少。此外, 多结点单元形函数阶次较高, 能够更好地模拟梁单元的变形形态。建议在工程实际中, 应根据计算精度 和效率的要求, 选择较少划分的多结点铁木辛柯梁单元或较多划分的2结点铁木辛柯梁单元。

　　 (2) 在弹性阶段, 采用静力凝聚可以很好地模拟单元性能, 并缩减自由度。而当结构进入强非线性时, 采用静力凝聚的结果不尽人意, 因此, 在采用多结点铁木辛柯纤维梁计算弹塑性动力时程分析时, 应慎重采用静力凝聚方法, 而采用内部结点出口的方式则结果良好。

　　 (3) 在截面板件较厚的梁单元进行纤维划分时, 若不考虑板厚影响, 会使计算截面精度下降, 造成结果失真, 因此在实际工程中需要注意。

　　 目前, 本文所述单元已在北京市建筑设计研究院有限公司自主研发的设计软件Paco中完整实现, 验证结果令人满意。