最新的《钢结构设计标准》 (GB 50017—2017) ^[1] (简称钢标) 针对结构分析和稳定性设计, 阐述了四种精度不一的分析和设计方法, 如表1所示。文献[2]也提出了类似的表格。这四种分析和设计方法主要从几何非线性和材料非线性的角度进行分类, 共考虑3个因素。其中, 整体二阶效应和构件二阶效应属于几何非线性的范畴, 而弹塑性效应属于材料非线性的范畴。

不同分析类型的非线性效应 表1

　　 注:“×”表示在分析或设计中, 该非线性效应不被考虑;“√”表示在分析或设计中, 该非线性效应被考虑。

　　 一阶弹性法在结构整体分析中不考虑几何非线性和材料非线性, 而是在分析结束之后在设计中通过采用有效长度或弯矩放大等方法来考虑非线性效应。该方法实质上将分析和设计割裂成两个独立的部分, 它一方面会造成分析过程中难以得到结构的真实响应, 另一方面因为不能准确得到计算长度从而为设计的安全性带来隐患。

　　 二阶弹性P-Δ或P-Δ-δ方法在分析阶段分别在不同层面上考虑了几何非线性所带来的弯矩放大效应, 从而减少了在设计阶段相关方面所采用的各种保守假设。特别是二阶弹性P-Δ-δ方法不仅在结构体系层面考虑了由于几何非线性产生的二阶效应, 而且从构件层面反映了由于轴力所带来的附加弯矩。如果在该方法的分析过程中引入几何初始缺陷, 就能在分析阶段直观地反映构件的稳定性情况。前3种方法 (一阶弹性、二阶弹性P-Δ或P-Δ-δ) 的共同特点都是在分析过程中不容许材料的弹塑性发展和内力重分布, 因此设计的时候应采用第一塑性铰原则。需要指出的是, 二阶弹性P-Δ-δ方法除了材料保持弹性, 其他要求都跟直接分析法的要求是一样的, 故为直接分析法的一种特例。

　　 当结构本身有较高的冗余度时, 可以考虑材料的非线性发展, 以获得更经济的设计。直接分析设计法在分析过程中显式地考虑材料的弹塑性效应, 会更加真实地模拟结构的工作状态。特别是在地震分析时, 采用弹塑性分析可以更加客观地反映结构在地震作用下的响应, 而构件的极限状态可根据设计目标及构件在整个结构中的作用来确定。直接分析设计法将结构分析和构件设计有机地结合在一起。直接分析设计法应考虑整体和构件二阶效应、结构和构件的初始缺陷、节点连接刚度^[3]、容许材料的弹塑性发展和内力重分布。从数值模拟的角度, 整体和构件二阶效应需要模拟构件的有限单元具备考虑局部二阶效应的能力^[4,5], 结构整体初始缺陷可通过偏移节点坐标或施加名义荷载来实现, 构件的初始缺陷需要使用二阶梁柱单元或多根梁柱单元来模拟构件。材料的非线性特性则应在有限元单元层面考虑。本文将重点探讨直接分析设计法中的二阶弹塑性分析。

　　 数值分析时, 可采用三维的实体单元、二维壳单元或一维梁柱单元进行模拟。其中, 实体单元和壳单元能够更好地模拟结构的弹塑性响应, 但建模复杂且计算时间长, 不利于工程应用, 因此这两类单元主要用于针对试验进行参数化分析的科学研究。梁柱单元则具有建模方便 (特别在采用一构件一单元的方法时) 、计算量小等特点, 与规范设计公式配套, 适合实际工程的分析设计。

　　 梁柱单元主要采用塑性铰法和塑性区法来进行二阶弹塑性分析。塑性铰法通过定义类似M-θ函数关系, 间接地描述构件的非线性弹塑性特性, 计算量小。但是它的计算精度相对偏低, 且这种类似M-θ的函数关系一般只是针对单轴受弯的情况, 而忽略了双轴耦合或轴力的影响。塑性区法则采用纤维截面显式地表征截面和构件的屈服情况, 具有更高的计算精度, 且能够反映P-M₂-M₃三者的耦合效应。对于钢构件, 残余应力会影响构件的刚度变化。纤维截面则能显式地模拟残余应力, 并具有高度的适用性, 能表征不同的残余应力分布模式。

　　 本文将针对直接分析设计法, 介绍一种基于柔度法的二阶梁柱单元。该单元采用力的平衡关系作为形函数, 具有较高的精度, 同时, 在单元层面考虑P-δ效应, 并在单元推导过程中直接使用函数表征构件的几何初始缺陷, 满足规范针对几何初始缺陷的相关要求。运用基于纤维截面技术的塑性区法, 能够显式地反映钢构件残余应力情况, 较为精细地模拟构件沿截面方向和沿构件方向的屈服过程。最后, 本文选用了几个经典算例验证了该单元的性能。

　　 根据假设的未知变量不同, 梁柱单元一般可基于变形假设推导而成, 简称刚度法, 或基于力的平衡假设推导而成, 简称柔度法。刚度法单元的形函数一般为多项式, 而多项式的次数决定了单元的精度, 如Chan 和Zhou ^[6]推导的PEP单元采用的是五次多项式, 具有很高的精度。对于刚度法单元的形函数, 其采用的多项式次数越多, 其公式的复杂性也越高。柔度法单元的刚度矩阵一般由数值积分而成, 其精度一般取决于积分点的数目, 而数目的增加并不影响公式的推导及其复杂性。

　　 下面将介绍文献[7]推导的一种基于柔度法的高性能梁柱单元 (FBMI) 。该单元能够直接反映构件的几何初始缺陷和材料初始缺陷, 并具有很高的精度, 可更加真实地表征构件的实际形态, 满足直接分析法适用于实际工程设计。

　　 以变形为变量的刚度法单元, 一般由最小势能原理推导而成, 而FBMI是由Hillinger-Reissner (HR) 变分方法进行推导。HR变分方法的表达式为:

　　 $\mathbf{\Pi }{}_{\text{Η}\text{R}}(\sigma ,u)={\displaystyle \underset{\Omega }{\int }\{}\epsilon (x,y,z)\sigma -\chi (\sigma )\}\text{d}\Omega +\mathbf{\Pi }{}_{\text{e}\text{x}\text{t}}(u)(1)$

　　 式中:u为位移;σ为应力;ε为应变;χ (σ) 为余能密度;Ω为构件体积;x, y和z为局部坐标系中的坐标。

　　 通过对式 (1) 变分可得:

　　 $\delta \mathbf{\Pi }{}_{\text{Η}\text{R}}(S,u)=\delta {}_{\text{S}}\mathbf{\Pi }{}_{\text{Η}\text{R}}+\delta {}_{\text{u}}\mathbf{\Pi }{}_{\text{Η}\text{R}}=0(2)$

　　 式中S为合应力场。

　　 取式 (2) 的弱形式, 可得力的平衡方程和变形协调方程:

　　 $\begin{array}{l}\delta {}_{\text{u}}\mathbf{\Pi }{}_{\text{Η}\text{R}}=0(3)\\ \delta {}_{\text{S}}\mathbf{\Pi }{}_{\text{Η}\text{R}}=0(4)\end{array}$

　　 力的平衡方程式 (3) 可以进一步表述为:

　　 $\delta {}_{\text{u}}\mathbf{\Pi }{}_{\text{Η}\text{R}}={\displaystyle \underset{L}{\int }\mathit{S}}{}^{\text{Τ}}\left\{\begin{array}{c}(\delta u^{\prime }+v^{\prime }\delta v^{\prime }+w^{\prime }\delta +\\ v^{\prime }{}_0\delta v´+w^{\prime }{}_0\delta w´)\\ \delta v^{″}\\ -\delta w^{″}\\ \delta {\psi }^{\prime }\end{array}\right\}\text{d}x-\overline{\mathit{{\rm P}}}{}^{\text{Τ}}\delta \mathit{D}=0(5)$

　　 式中:v₀, w₀为构件初始缺陷;v, w为侧向位移;ψ为扭转角度;P为构件端部力;D为节点位移;L为构件长度。

　　 通过带入边界条件, 可得截面内力与杆端力之间的关系为:

　　 $\mathit{S}(x)=\left\{\begin{array}{c}{\rm N}(x)\\ {\rm M}{}_{\text{z}}(x)\\ {\rm M}{}_{\text{y}}(x)\\ {\rm T}(x)\end{array}\right\}=\mathit{b}(x)\mathit{{\rm P}}(6)$

　　 式中:N为轴力;M_z, M_y为弯矩;T为扭转力;b为系数矩阵。

　　 变形协调方程式 (4) 可以展开为:

　　 $\begin{array}{l}\delta {}_{\text{S}}\mathbf{\Pi }{}_{\text{Η}\text{R}}=\\ {\displaystyle \underset{L}{\int }\delta }\mathit{S}{}^{\text{Τ}}\left[\left\{\begin{array}{c}(u´+\frac{1}{2}v{}^{2^{\prime }}+\frac{1}{2}w{}^{2^{\prime }}+\\ v^{\prime }v{}_{0^{\prime }}+w^{\prime }w{}_{0)}´\\ v^{″}\\ -w^{″}\\ {\psi }^{\prime }\end{array}\right\}-\frac{\partial \chi (\mathit{S})}{\partial \mathit{S}}\right]\text{d}x=0(7)\end{array}$

　　 根据虚功原理可得:

　　 ${\displaystyle \underset{L}{\int }\delta }\mathit{S}(x){}^{\text{Τ}}\mathit{d}(x)\text{d}x=\delta \mathit{{\rm P}}{}^{\text{Τ}}\mathit{D}(8)$

　　 最后, 可得构件端部位移与截面变形之间的协调关系为:

　　 $\mathit{D}={\displaystyle \underset{L}{\int }\mathit{b}{}^{*}(x)}{}^{\text{Τ}}\mathit{d}(x)\text{d}x(9)$

　　 式中:d为杆件内部变形;b^*为系数矩阵。

　　 联合平衡方程和协调方程, 可得单元柔度矩阵为:

　　 $\begin{array}{c}\mathit{F}{}_{\text{e}}=\frac{\partial \mathit{D}}{\partial \mathit{{\rm P}}}\hfill \\ ={\displaystyle {\int }_L\{\mathit{b}{}^{*}(x)}{}^{\text{Τ}}\mathit{f}{}_{\text{s}}(x)[\mathit{b}(x)+\mathit{h}(x)]+\mathit{g}(x)\}\text{d}x\hfill \end{array}(10)$

　　 式中f_s为截面柔度矩阵。其中, h (x) 与g (x) 为:

　　 $\begin{array}{l}\mathit{h}(x)={\rm P}{}_1\left[\begin{array}{c}0\\ \mathit{V}(x)\\ -\mathit{W}(x)\\ 0\end{array}\right](11)\\ \mathit{g}(x)=\frac{1}{2}\kappa {}_{\text{z}}\left[\begin{array}{c}\mathit{V}(x)\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]-\frac{1}{2}\kappa {}_{\text{y}}\left[\begin{array}{c}\mathit{W}(x)\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right](12)\end{array}$

　　 其中, κ_z, κ_y为曲率;V (x) 和W (x) 为:

　　 $\begin{array}{l}\mathit{V}(x)=\frac{\partial v(x)}{\partial \mathit{{\rm P}}}(13)\\ \mathit{W}(x)=\frac{\partial w(x)}{\partial \mathit{{\rm P}}}(14)\end{array}$

　　 式中v (x) , w (x) 包含有几何初始缺陷v₀, w₀。

　　 在1.1节单元构件的推导过程, 显式地包含了构件的几何初始缺陷v₀, w₀。

　　 在钢标中, 构件几何初始缺陷的模式可以定义为:

　　 $\delta {}_0=e{}_0\sin \frac{\text{π}x}{l}(15)$

　　 式中:δ₀为离构件端部x处的初始变形量, 也即v₀和w₀;e₀为构件中点处的初始变形量;x为距离构件端部的距离;l为构件的总长度。

　　 图1显式地描述构件的初始缺陷。

　　 钢标同时提供了另外一种方法模拟初始缺陷——假想均布荷载, 如图2所示。假想在分析计算之前, 根据初始缺陷的方向和大小, 将均布荷载施加于构件。与在单元层面显式表征初始缺陷相比, 假想荷载方法将增加编写计算程序的复杂性, 并增加计算量。而1.1节推导的单元显式地考虑初始缺陷, 能满足钢标的相关要求。

　　 在直接分析法中, 如果采用塑性区法考虑材料非线性, 构件的几何初始缺陷e₀应按不小于l/1 000的出厂加工精度取值, 并考虑初始残余应力。

　　 材料初始缺陷一般指的是广泛存在于热轧或焊接钢构件中的残余应力。残余应力主要是因为轧制、焊接以及冷成型等工艺造成的。残余应力在构件中表现为应力自平衡, 应力分布模式与截面形状具有极大的相关性。例如, 工字形截面或H形截面的翼缘中间部位表现为拉应力, 而其他部位表现为压应力;方形截面的四个角部表现为拉应力, 而其他部位表现为压应力。图3为两种常用截面类型的残余应力分布模型。

　　 残余应力作为一种初始的力学缺陷, 会造成构件提前屈服以及裂纹的出现, 也是影响钢结构轴压构件稳定性能的因素之一。

　　 钢标提供两种方法来考虑残余应力的对构件承载力的影响。第一种方法是使用构件的综合缺陷代表值, 同时表征几何初始缺陷和残余应力。当采用直接分析而不考虑材料弹塑性发展时, 构件综合缺陷代表值如表2所示。第二种方法利用塑性区法, 考虑初始残余应力。本文介绍的柔度法单元, 采用纤维截面法能显式地在单元层次上考虑残余应力。

　　 在直接分析法中, 如果不考虑材料弹塑性发展时, 结构分析应限于第一个塑性铰的形成, 且不允许内力重分布;如果考虑材料弹塑性发展, 直接分析法宜采用塑性铰法或塑性区法。钢标推荐弹塑性分析中材料的本构关系采用理想弹塑性模型, 屈服强度可取强度设计值, 弹性模量应取标准值。对于一些特殊分析, 应采用相应的本构模型, 如连续倒塌分析应采用考虑应变率影响的应力-应变关系;抗火分析应考虑结构材料在高温下的应力-应变关系。钢标针对材料弹塑性发展, 在构造上做了相应的限定, 如塑性铰形成的区域, 构件和节点应有足够的延性保证以便内力重分布, 容许一个或多个塑性铰产生;钢结构构件截面应为双轴对称截面或单轴对称截面, 塑性铰处截面宽厚板等级应为S1, S2级, 其出现的界面或区域应保证有足够的转动能力。

　　 塑性铰法通常在基于刚度法的梁柱单元中使用, 相当于在梁柱单元的两端附着两个零长度的弹簧, 从而组成的一个超级单元, 如图4所示。该做法与考虑节点半刚性时的梁柱做法类似, 不同的是, 在构件端部的内力未使构件端部出现塑性时, 该弹簧的刚度为无穷大;而当构件端部逐渐由弹性进入全塑性, 弹簧的刚度由无穷大逐步变为0。

　　 构件端部逐渐由弹性进入全塑性的变化过程一般由假设的函数模型所表示。常用的塑性铰模型如图5所示。

　　 塑性铰模型存在两个方面的问题。一方面, 从图5中可知, 塑性铰模型过于简化, 一般不能较为真实地反映构件进入塑性的过程。另一方面, 塑性铰模型一般只考虑绕单轴的弯矩效应, 而忽略了轴力以及绕另一个轴的弯矩的耦合效应, 从而限制了塑性铰的适用范围。

　　 塑性区法一般被基于柔度法的梁柱单元所采用, 该方法使用纤维来离散横截面, 使用数值积分点沿构件长度方向来离散构件, 如图6所示。沿构件方向的常用的数值积分方法有:Gauss-Lobatto 法、Gauss-Legendrei 法和Newton-Cotes integration法。文献[8]详细讨论了各种数值积分方法的精度。总体而言, 基于柔度法的梁柱单元采用数值积分, 当积分点达到7个就可以得到较高的计算精度, 能满足实际工程设计需要。

　　 图6的横截面由纤维截面表征, 而纤维截面中每根纤维的应力均由单独的参数存储, 具有独立性。对于钢材而言, 常用的材料非线性本构模型有:理想弹塑性模型和Giuffré-Menegotto-Pinto模型, 图7为钢材非线性本构模型示意图。理想的弹塑性模型相对Giuffré-Menegotto-Pinto模型, 具有简单、效率高的特点, 在工程设计和理论研究中都较为常用。Giuffré-Menegotto-Pinto模型有较为平滑的弹性-塑性过渡段, 有利于非线性的收敛, 一般常用于理论研究中。

　　 直接分析设计法对于弹塑性效应的考虑, 主要是通过在非线性分析中考虑材料的弹塑性特性实现。以塑性区法为例, 图8 (虚线框部分) 展示了材料非线性在有限元单元迭代计算中所处的位置。

　　 首先, 根据材料类型选用适当的非线性的材料本构关系 (图7) 。纤维截面中每根纤维使用独立的本构模型 (图6) , 从而纤维截面模型能较为光滑地表征截面从弹性到完全塑性的过程。在为每根纤维赋予材料模型之后, 可以根据残余应力分布情况, 为每根纤维设置初始应力, 即显式表征残余应力。

　　 然后, 截面的合力和刚度都通过沿纤维截面的积分所得。在沿截面积分过程中, 材料的弹塑性效应在两方面体现:1) 每根纤维均采用非线性本构模型, 并各自有独立的数据记录每根纤维的状态;2) 由沿纤维截面的数值积分所体现的非线性。随着非线性迭代以及荷载的加大, 纤维截面的状态会不断被更新, 有些纤维逐渐进入塑性, 造成截面刚度K_s的弱化, 进而减弱了该构件及整个结构体系的刚度。

　　 图9为两层钢框架结构, 该算例具有较强的几何非线性和材料非线性特点, 最早被Ziemian等 ^[9]所研究。Du等^[7]使用该算例验证了其推导的文中所介绍的基于柔度法的梁柱单元, 本文将从弹塑性分析的角度解析该算例。图9中体现了框架的平面图和加载情况, 荷载组合为1.2倍的恒荷载+1.6倍的活荷载。所有的构件均采用相同的钢材, 弹性模量为199 810 MPa, 屈服应力为248 MPa。1层柱底的边界条件为铰接, 所有的梁柱节点假设为刚性连接。采用7种模型模拟该算例, 见表3。

　　 表3中, M1～M3采用的初始缺陷模式为结构一阶屈曲模态, 如图10所示。所有的材料模型均采用理想的弹塑性模型。所有模型在计算中, 记录图9右上角柱顶 (圆圈处) 侧向位移情况。

　　 图11为7种模型的荷载系数-位移的关系曲线。从图11可知, M4与M5采用塑性铰法, 当荷载系数达到1.0之后未能出现位移方向的反转, 而是呈现出发散的态势。出现该情况的原因可能是由于塑性铰模型的不连续, 不能从弹性状态平滑地过渡到全塑性状态, 从而引起了数值上的不稳定。而M1～M3, M6及M7都采用了塑性区法, 则统一得到了相对平滑且稳定的模拟结果。M2和M3分别从不同层面考虑几何初始缺陷, 很好地模拟了初始缺陷所带来的刚度弱化效应。

　　 图12为M3计算得到的内力图以及构件屈服情况。从该图可知, 部分梁柱节点处的截面有90%以上的纤维处于屈服状态。

　　 直接分析设计法通过全面考虑二阶弹性P-Δ和P-Δ-δ效应、结构体系和构件初始缺陷、节点连接刚度和材料非线性, 获取各种荷载设计值下的内力和标准值下的位移, 直接反映构件真实响应, 从而不需要计算长度法进行构件受压稳定承载力验算。针对材料的弹塑性发展和结构的内力重分布, 需要采用二阶弹塑性分析。本文介绍了一种基于柔度法的高性能梁柱单元, 该单元能够满足钢标对初始缺陷和几何非线性要求的同时, 采用塑性区法以及纤维截面技术来模拟材料非线性, 满足钢标中直接分析法中的二阶弹塑性分析。得到的主要结论为:

　　 (1) 材料非线性在直接分析法中占有重要的地位, 是将直接分析法从适用于静力设计推广到动力时程分析中的一个关键因素。

　　 (2) 与工程实际中普遍使用的塑性铰法相比, 塑性区法能够更好地反映结构的弹塑性响应, 但该方法也存在计算量大的问题。

　　 (3) 一种基于柔度法的二阶梁柱单元在几何非线性和材料非线性两方面都具有优异的表现, 并能满足钢标中直接分析设计法的相关要求, 能够在实际的工程分析设计中使用。