变截面构件材料利用率高、自重轻, 可满足现代建筑外观及使用功能的要求, 被大量用于各种结构形式中^[1,2,3,4]。此外, 自动化焊接与建筑信息模型 (BIM) 技术的发展与普及大大降低了构件的制造及施工成本, 推动了变截面构件在结构工程中的应用。然而, 现行设计方法存在诸多问题, 给结构安全带来了一定隐患。

　　 目前, 变截面构件的内力和变形计算以一阶线性分析方法为主。由于现有梁柱单元主要针对等截面构件, 工程中通常采用两种方法模拟变截面单元。第一种方法为阶梯单元法, 方法采用一系列等截面梁柱单元来代替一根变截面构件。根据经验, 为准确反映构件的受力与变形, 一根变截面构件需至少划分为20根等截面单元^[5,6]。该方法可以较精确地模拟变截面构件的力学行为, 但建模和计算工作量大, 不适用于结构形式较复杂的工程。第二种方法为近似刚度法, 方法假定构件弯曲刚度沿长度近似为线性、二次或三次变化函数^[7]。该方法在建模与分析过程中与等截面构件并无区别, 建模及计算效率高, 对截面形式简单的构件, 如实心矩形, 计算结果较准确, 但不适用于复杂截面形式^[8]。

　　 传统设计方法通过基于试验结果的经验公式评估构件的稳定性。设计公式往往基于有效长度法, 采用等效惯性矩或等效长度对等截面构件设计公式进行修正, 如美国变截面构件设计指南^[9]和中国《门式刚架轻型房屋结构技术规范》 (GB 51022—2015) ^[10]。规范中的设计公式, 仅仅适用于变高度构件的设计。对边界条件、变截面形式与试验条件相差较大的构件, 传统设计方法无法正确计算构件的承载力。

　　 因此, 本研究将直接分析法理论应用于变截面构件的设计。直接分析法为基于有限元理论的数值方法, 该方法在结构数值分析过程中直接考虑各种与结构稳定有关因素的影响, 因此无需采用传统的有效长度假定^[11,12,13], 尤其适用于形式复杂的结构。目前, 由于一单元一构件的模拟方式在实际工程设计中优势明显, 已成为直接分析法中的首选建模方式^[14]。与传统梁柱单元不同^[15,16,17], 本文针对一单元一构件的建模方式, 推导出包含最不利几何初始缺陷影响的高阶变截面梁柱单元, 提出变截面构件承载力计算公式。验证算例表明, 提出的方法可靠准确, 可用于各种变截面形式的分析与设计。

　　 构件的几何缺陷会增大杆件P-δ效应, 是产生非线性行为的重要因素。工程中通常采用简支柱在轴压力下的一阶屈曲模态作为构件的最不利初始形状。因此, 等截面构件的初始弯曲可用二次抛物线函数描述。然而, 对于变截面构件, 该函数无法代表其最不利的初始缺陷 (图1) 。

　　 因此, 本研究提出用于描述变截面构件初始缺陷y_m0 (ζ) 的多项式如下:

　　 $y{}_{\text{m}0}(\zeta )=v{}_{\text{m}0}{\displaystyle \sum _{i=0}^{5}C}{}_i\zeta {}^i(1)$

　　 式中:v_m0为初始缺陷幅值;$\zeta =\frac{x}{L}$, x为截面在x坐标轴中的位置, L为构件长度, 根据规范^[18]规定选取;C_i为常系数, 其值见表1。

　　 为简化推导过程和单元表达式, 采用以下假定:1) 大变形小应变假设;2) 忽略剪切和翘曲变形的影响;3) 不考虑局部屈曲;4) 节点力均为保守力;5) 构件高度及宽度沿长度方向均为线性变化, 而翼缘及腹板厚度保持不变。

　　 局部坐标系下, 单元节点力和位移如图2所示。y轴和z轴分别表示截面的弱轴和强轴。单元由一个可自由移动的内部节点和两个子单元组成, 并包含初始缺陷的影响。在结构分析过程中, 该内部节点可以自由移动至最危险截面位置, 以捕捉构件最薄弱处的内力与变形。单元的变形可以通过节点自由度表示, 以矩阵形式可表示为:

　　 $\widehat{\mathit{U}}=\mathit{{\rm N}}(\zeta )\widehat{\mathit{u}}(2)$

　　 其中:

　　 $\begin{array}{l}\widehat{\mathit{U}}{}^{\text{Τ}}=\left[u{}_{\text{a}},v{}_{\text{y}\text{L}},v{}_{\text{y}\text{R}},u{}_{\text{t}},v{}_{\text{z}\text{L}},v{}_{\text{z}\text{R}}\right](3)\end{array}$

$\begin{array}{l}\\ \widehat{\mathit{u}}{}^{\text{Τ}}=[e,\theta {}_{\text{y}11},\theta {}_{\text{z}11},\theta {}_{\text{t}},\theta {}_{\text{y}12},\theta {}_{\text{z}12},\delta {}_{\text{z}},\delta {}_{\text{y}},\theta {}_{\text{y}21},\theta {}_{\text{z}21},\theta {}_{\text{y}22},\theta {}_{\text{z}22}](4)\end{array}$

$\begin{array}{l}\\ \begin{array}{c}\mathit{{\rm N}}(\zeta )=[\begin{array}{cccccccccccc}{\rm N}{}_1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\end{array}\end{array}$

$\begin{array}{l}\begin{array}{c}\begin{array}{c}\mathrm{}\\ 0& 0\end{array}\end{array}\end{array}$Ν200Ν30Ν40000

0Ν200Ν30Ν400000

000Ν100000000

0000000Ν50Ν60Ν7

000000Ν50Ν60Ν70]          (5)

　　 式中:$\widehat{\mathit{U}}$为单元变形矩阵;$\widehat{\mathit{u}}$为节点自由度矩阵;u_t为单元总扭转位移;δ_z, δ_y分别为节点位移;u_a为单元轴向位移函数;v_yL 和v_zL分别为左侧子单元沿y轴和z轴的侧向位移函数; v_yR 和v_zR分别为右侧子单元沿y轴和z轴的侧向位移函数; θ_y11, θ_z11, θ_y12, θ_z12, θ_y21, θ_z21, θ_y22, θ_z22为节点转角, 其意义见图2;e为单元轴向总位移;θ_t为单元总扭转角;N (ζ) 为形函数矩阵, 其各元素具体表达式为:

　　 $\begin{array}{l}{\rm N}{}_1=\frac{1}{2}+\zeta (6)\end{array}$

$\begin{array}{l}\\ {\rm N}{}_2=\frac{2L(1+2\zeta )(\zeta -\xi ){}^2}{(1+2\xi ){}^2}(7)\end{array}$

$\begin{array}{l}\\ {\rm N}{}_3=\frac{L(1+2\zeta ){}^2(\zeta -\xi )}{(1+2\xi ){}^2}(8)\end{array}$

$\begin{array}{l}\\ {\rm N}{}_4=\frac{(1+2\zeta ){}^2(1-4\zeta +6\xi )}{(1+2\xi ){}^3}(9)\end{array}$

$\begin{array}{l}\\ {\rm N}{}_5=-\frac{2L(1-2\zeta )(\zeta -\xi ){}^2}{(1-2\xi ){}^2}(10)\end{array}$

$\begin{array}{l}\\ {\rm N}{}_6=\frac{L(1-2\zeta ){}^2(\zeta -\xi )}{(1-2\xi ){}^2}(11)\end{array}$

$\begin{array}{l}\\ {\rm N}{}_7=\frac{(1-2\zeta ){}^2(1-4\zeta +6\xi )}{(2\xi -1){}^3}(12)\end{array}$

$\begin{array}{l}\\ -\frac{1}{2}<\xi <\frac{1}{2}(13)\end{array}$

　　 在单元矩阵推导过程中, 截面参数沿长度方向的变化采用精确表达式。构件的截面高度B (x) 及宽度D (x) 可表示为:

　　 $\begin{array}{l}B(x)=\frac{B{}_{\text{R}}-B{}_{\text{L}}}{L}\left(x+\frac{L}{2}\right)+B{}_{\text{L}}(14)\end{array}$

$\begin{array}{l}\\ D(x)=\frac{D{}_{\text{R}}-D{}_{\text{L}}}{L}\left(x+\frac{L}{2}\right)+D{}_{\text{L}}(15)\end{array}$

　　 式中:B_L, D_L分别为构件左端截面宽度和高度;B_R, D_R分别为构件右端截面宽度和高度。

　　 工字形截面惯性矩I_y (x) , I_z (x) 和横截面面积A (x) 沿杆轴向变化的表达式为:

　　 $\begin{array}{l}{\rm I}{}_{\text{y}}(x)=\frac{2{\rm T}{}_{\text{f}}B(x){}^3+\left[D(x)-2{\rm T}{}_{\text{f}}\right]t{}_{\text{w}}{}^3}{12}(16)\end{array}$

$\begin{array}{l}\\ {\rm I}{}_{\text{z}}(x)=\frac{t{}_{\text{w}}\left[D(x)-2{\rm T}{}_{\text{f}}\right]{}^3}{12}+\frac{B(x){\rm T}{}_{\text{f}}{}^3}{6}+\\ 2{\rm T}{}_{\text{f}}B(x)\left[\frac{D(x)}{2}-\frac{{\rm T}{}_{\text{f}}}{2}\right]{}^2(17)\end{array}$

$\begin{array}{l}\\ A(x)=2B(x){\rm T}{}_{\text{f}}+\left[D(x)-2{\rm T}{}_{\text{f}}\right]t{}_{\text{w}}(18)\end{array}$

　　 式中:I_y (x) 和I_z (x) 分别表示绕弱轴和强轴惯性矩;T_f和t_w分别为截面的翼缘和腹板厚度。

　　 单元刚度矩阵采用势能驻值原理进行推导。单元的总势能Π可以表示为:

　　 $\Pi =U-V(19)$

　　 式中:U为应变能;V为节点力做功。

　　 考虑初始缺陷的影响, 空间三维高阶变截面梁柱单元的应变能的表达式为:

　　 $\begin{array}{l}U=\frac{1}{2}{\displaystyle {\int }_{-L/2}^{L/2}E}A(x)\dot{u}{}^2\text{d}x+\frac{1}{2}{\displaystyle {\int }_{-L/2}^{L/2}E}{\rm I}{}_{\text{y}}(x)\stackrel{\cdot \cdot }{{\displaystyle v}}{}_{\text{y}}{}^2\text{d}x+\\ \frac{1}{2}{\displaystyle {\int }_{-L/2}^{L/2}E}{\rm I}{}_{\text{z}}(x)\stackrel{\cdot \cdot }{{\displaystyle v}}{}_{\text{z}}{}^2\text{d}x+\frac{1}{2}{\displaystyle {\int }_{-L/2}^{L/2}G}\stackrel{—}{{\displaystyle J}}\dot{u}{}_{\text{t}}^2\text{d}x+\\ \frac{1}{2}{\displaystyle {\int }_{-L/2}^{L/2}{\rm P}}(\dot{v}{}_{\text{y}}{}^2+2\dot{v}{}_{\text{y}}\dot{y}{}_{\text{m}0\text{y}})\text{d}x+\\ \frac{1}{2}{\displaystyle {\int }_{-L/2}^{L/2}{\rm P}}(\dot{v}{}_{\text{z}}{}^2+2\dot{v}{}_{\text{z}}\dot{y}{}_{\text{m}0\text{z}})\text{d}x(20)\end{array}$

　　 式中:E和G分别为杨氏及剪切模量;P为轴力;$\stackrel{—}{{\displaystyle J}}$为截面等效扭转常数;$\dot{v}$和$\stackrel{\cdot \cdot }{{\displaystyle v}}$为函数v的一阶和二阶导数;下角标y和z分别代表对应的截面形心主轴;$\dot{u}$和$\dot{u}{}_t$分别为轴向位移函数u和扭转位移函数u_t的一阶导数。

　　 外力做功:

　　 $V={\displaystyle \sum _{i=1}^{12}F}{}_{i}u{}_i(21)$

　　 式中:F_i和u_i分别为节点力和位移;i为对应的自由度编号。

　　 根据势能驻值原理, 式 (19) 的一阶导数为单元平衡方程:

　　 $\frac{\delta \Pi }{\delta u{}_i}=\frac{\partial U}{\partial u{}_i}+F{}_i=0(i=1\sim 12)(22)$

　　 式 (22) 为单元的割线刚度关系, 即节点力表达式。

　　 对单元平衡方程求导, 即可得单元的刚度:

　　 $k{}_{ij}=\frac{\partial F{}_i}{\partial u{}_j}=\frac{\partial U}{\partial u{}_i\partial u{}_j}(i,j=1\sim 12)(23)$

　　 式中i, j为对应自由度编号。

　　 单元刚度矩阵由3部分组成:

　　 $\mathit{{\rm K}}{}_{\text{E}}=\mathit{{\rm K}}{}_{\text{L}}+\mathit{{\rm K}}{}_{\text{G}}+\mathit{{\rm K}}{}_{\text{S}}(24)$

　　 式中K_E, K_L, K_G和K_S分别为单元总刚度矩阵、线性刚度矩阵、几何刚度矩阵和内部节点转动刚度。

　　 节点力和刚度矩阵的具体表达式见文献[2]。

　　 与等截面构件类似, 直接分析法中, 变截面构件轴压柱的稳定性ψ可通过下式计算:

　　 $\psi =\max \left[\frac{{\rm P}}{{\rm P}{}_{\text{c}}(x)}+\frac{{\rm P}v{}_{\text{y}}(x)}{{\rm M}{}_{\text{c}\text{z}}(x)}+\frac{{\rm P}v{}_{\text{z}}(x)}{{\rm M}{}_{\text{c}\text{y}}(x)}\right](25)$

　　 式中:max[F (x) ]表示函数F (x) 的最大值; v (x) 为包含初始缺陷影响的构件侧向位移;P_c (x) =fA (x) , M_c (x) =f W (x) 分别为截面轴力和弯矩承载力, f为材料设计强度, W (x) 为截面模量;下角标y和z分别代表对应的截面形心主轴。

　　 式 (25) 为构件内力、变形和截面属性的函数。

　　 当ψ值不大于1.0时, 可认为构件在该荷载条件下是稳定的, 否则应对构件尺寸重新设计。

　　 为验证提出的缺陷模型的准确性, 将由式 (1) 计算出的曲线与实体单元计算出的一阶屈曲模态进行对比。对于不同变截面形式的工字形变截面柱, 对比结果如图3, 4所示, 图中ζ表示截面位置, y_m0/ν_m0表示侧向位移比。

　　 由以上对比结果可知, 所提出的缺陷模型可较精确地模拟构件在不同变截面形式下的最不利初始构形。因此, 式 (1) 可应用于一单元一构件的变截面构件直接分析法。

　　 单元精度决定了分析结果的可靠度, 优秀的分析单元应具有分析精度高、计算速度快的特点。算例采用一单元一构件的模拟方式, 并与传统阶梯单元法进行对比, 通过对简支变截面构件在不同荷载条件下的分析, 验证单元的精度和效率。构件在小端和大端的截面尺寸分别为500mm×1 000×30mm×25mm和1 000mm×500mm×30mm×25mm, 构件长度为20m。弹性模量为205GPa, 泊松比为0.3。荷载施加方式、荷载数值及构件轴向变形计算结果见图5。计算结果汇总与对比见表2。

　　 由图5可知, 构件在所施加荷载作用下表现出明显的非线性行为, 其中, 变截面高阶单元的最大计算误差为3.92%, 平均误差为1.74% (表2) , 说明该一单元一构件方法能够精确模拟变截面构件的非线性行为。与传统阶梯单元方法相比, 该方法计算效率高, 建模简便, 可在实际工程中广泛应用。

　　 直接分析法的核心是在结构分析计算过程中直接考虑对稳定有影响的各个因素, 例如构件和结构的初始缺陷。为验证高阶变截面梁柱单元在直接分析法中的精度, 本算例研究了包含初始缺陷的轴压柱在不同变截面形式及边界条件下的变形行为, 并与传统有限元方法进行了对比。构件截面尺寸及初始缺陷见图6, 柱高均为5m。传统方法通过施加一阶屈曲模态, 对构件的初始缺陷进行模拟。由两种方法计算出的构件轴向位移如图6所示。

　　 由图6可知, 该一单元一构件的模拟方法可准确考虑最不利初始缺陷对变截面构件的影响, 计算结果精度高, 可大大简化工程设计中的建模过程, 提高计算效率。

　　 为验证构件承载力计算公式 (式 (25) ) 的可靠性, 本算例根据美国钢结构规范ANSI/AISC 360-16^[12]的规定, 采用两种模拟方式, 即壳单元与高阶变截面梁柱单元, 对变截面构件弯曲失稳的承载力进行了对比计算。构件初始缺陷为L/1 000, 残余应力模型选用里海模型^[19]。

　　 为检验构件初始缺陷模拟方法的准确性, 将传统有限单元法及本文提出的直接分析法计算出的等截面柱失稳曲线与规范曲线进行了对比, 如图7所示。三条曲线基本一致, 可认为初始缺陷的模拟方法正确, 可应用于变截面构件的研究。

　　 对于不同变化程度的变截面构件, 以壳单元计算结果作为基准, 与直接分析法计算出的柱的失稳曲线对比见图8。

　　 由图8可知, 两种方法的计算结果相当, 表明该直接分析法计算结果可靠, 可应用于实际工程设计。

　　 将直接分析法理论推广应用至变截面构件, 推导了适用于一单元一构件模拟方式的欧拉-伯努利高阶变截面梁柱单元, 提出了变截面构件截面承载力验算公式。验证算例表明, 本文提出的变截面构件直接分析法计算结果准确, 具有建模简便, 计算效率高的特点, 可用于工程中各种复杂变截面形式的设计。