空间梁单元切线刚度矩阵的精确分析方法

刘树堂

摘要：为有效进行空间刚架结构后屈曲分析，提出一种新的空间梁单元切线刚度矩阵的精确分析方法。首先用直接法建立梁单元杆端力与杆端位移增量的关系式，然后根据矩阵微分理论求单元杆端力关于杆端位移的导数，最后在求导结果表达式中令杆端位移增量为0，即可得到梁单元切线刚度矩阵。对六层和二十层空间刚架结构进行了后屈曲分析。结果表明：所得的空间梁单元切线刚度矩阵具有足够精度，可有效用于大型空间刚架结构的后屈曲分析。

关键词：空间梁单元；切线刚度矩阵；空间刚架；非线性屈曲；弹塑性屈曲

中图分类号：TU393.3 文献标志码：A

0 引 言

在空间刚架结构的屈曲分析中，单元切线刚度矩阵的预测精度对分析结果的正确性有着关键性影响，特别是对于结构的后屈曲分析。若单元切线刚度矩阵预测精度不够，则很难通过荷载-位移平衡路径的极限点，因此也很难准确地确定结构的极限荷载。建立空间梁单元切线刚度矩阵有多种途径，如势能变分方法[1-2]、Timosheko梁柱理论[3]、弹性刚度矩阵与几何刚度矩阵的叠加方法[4]以及一些其他方法[5-7]。

势能变分方法[1]建立切线刚度矩阵常常要忽略一些高阶项，导致切线刚度矩阵的预测精度受到影响。Timosheko梁柱理论直接由平衡方程来导出切线刚度矩阵，平衡方程中的力与位移关系可用超越函数表示，分析方法简单，切线刚度矩阵更精确。

Yang等[4]首先采用虚功方法来建立梁单元的几何刚度矩阵，然后由弹性刚度矩阵与几何刚度矩阵的叠加建立切线刚度矩阵。文献[5]中基于Prandtl-Renss增量理论，推导出三维梁单元的弹塑性切线刚度矩阵。文献[6]中基于非线性问题的一般平衡方程和空间梁单元的非线性几何方程，推导应力-应变一般线弹性关系下的空间梁单元显式切线刚度矩阵，该刚度矩阵中包含了由初应力和初应变产生的初应力刚度矩阵。文献[7]中使用复合变量求导方法对内力近似求导得到切线刚度矩阵。上述这些方法所建立的空间梁单元切线刚度矩阵并没有通过大型算例来验证其有效性。

为有效进行空间刚架结构后屈曲分析，本文中提出一种新的空间梁单元切线刚度矩阵的分析方法。首先在变形状态下用直接法建立梁单元杆端力与杆端位移增量的关系式，然后根据矩阵微分理论求单元杆端力关于杆端位移增量的导数，最后令杆端位移增量为0，即可得到梁单元切线刚度矩阵。本文中笔者建立切线刚度矩阵时，没有忽略任何高阶项，因此所建立的切线刚度矩阵是精确的。通过六层和二十层空间刚架结构的后屈曲分析，验证了本文空间梁单元切线刚度矩阵的有效性。

1 空间梁单元杆端向量变换矩阵

空间刚架结构的精确非线性分析，需要利用节点空间大转动分析理论[8-10]。节点空间大转动情况的转动次序不满足交换律，即节点总量转角不能由转角增量线性叠加累积计算。

节点空间大转动分析需要建立3种局部坐标系：节点坐标系、单元坐标系和单元端面坐标系。这3种局部坐标系对整体坐标系的旋转矩阵分别称为节点方向矩阵、单元方向矩阵和单元端面方向矩阵。图1中给出了结构变形构形下空间梁单元的单元坐标系、单元端面坐标系和结构坐标系。

1.1 节点方向矩阵

节点坐标系是一个刚性固定于节点上的笛卡儿坐标系，节点发生转动时，该坐标系方向（节点方向）随之变化。通常令初始构形下的节点坐标系与整体坐标系平行，故初始节点方向矩阵q0为3阶单位矩阵，q0=I3，I3为3阶单位矩阵。节点发生转动后，节点方向矩阵q按照节点空间大转动的坐标旋转变换矩阵进行更新，即

式中：θi为转角矩阵θ的元素，i=1，2，3；θx，θy，θz分别为x，y，z方向的转角。

1.2 单元方向矩阵

单元坐标系以Oxeyeze表示（图1），单元坐标系的xe轴通过单元j端和k端截面形心。在初始构形下，单元坐标系的ye，ze轴分别与单元截面2个主轴方向平行。单元方向矩阵由单元坐标系的xe，ye，ze轴对整体坐标系x，y，z轴的方向余弦向量构成。初始构形下的单元方向矩阵r0表示为

式中：r01，r02，r03分别为初始构形下的单元坐标系xe，ye，ze轴对整体坐标系x，y，z轴的方向余弦向量。

在变形构形下，单元截面形心轴发生弯曲，单元2个端截面也发生了转动和位移，单元坐标系的方向发生了变化。在变形构形下，单元坐标系的xe轴对整体坐标系x，y，z轴的方向余弦向量r1为

式中：u为整体坐标系下单元j端和k端节点线位移集合向量，u=（ujx，ujy，ujz，ukx，uky，ukz）T；x为增量步开始时单元j端和k端节点坐标的集合向量；t1=（-I3，I3）。

变形构形下的单元方向矩阵r可表示为

式中：r1，r2，r3分别为变形构形下单元坐标系xe，ye，ze轴对整体坐标系x，y，z轴的方向余弦向量。

在变形构形下，单元坐标系ye，ze轴对整体坐标系x，y，z轴的方向余弦向量r2，r3需要由单元端面转动情况来确定。

1.3 单元端面方向矩阵

单元端面坐标系是一个固定于单元端截面上的笛卡儿坐标系（图1中的Oxjyjzj和Oxkykzk），用来定义单元端面的方向。

在初始构形下，单元端面坐标系与单元坐标系一致，单元端面方向矩阵即为初始构形下的单元方向矩阵r0。在变形构形下，节点发生转动，单元端面坐标系的方向发生改变，单元端面方向矩阵随之变化。由于单元端面与节点刚性连接，节点发生转动时，单元端面坐标系与节点坐标系同步转动，它们之间保持着固定变换关系。

设在上一增量步迭代结束时，单元j端和k端的节点在整体坐标系下的总量转角分别为θ0j和θ0k，此时单元j端和k端的单元端面方向矩阵分别为R（θ0j）r0和R（θ0k）r0。endprint

设在当前增量步上，单元j端和k端的节点分别发生整体坐标系下的转角增量为θj和θk。

1.4 确定单元方向矩阵的第2列和第3列

节点发生转动后，单元j端和k端的端面坐标系xj，xk轴的方向向量pj1，pk1分别为 将单元j端和k端的端面坐标系分别绕nj轴和nk轴转动γjnj和γknk角，使pj1和pk1与单元坐标系的xe轴重合。此时单元j端和k端的单元端面方向矩阵ej，ek分别变为

2 空间梁单元切线刚度矩阵

2.1 单元基本内力与单元基本变形的关系

单元基本内力增量向量与单元基本变形增量向量（均在单元坐标系下度量）的关系为

式中：V为单元基本变形增量向量，V=（v，，αyj，αyk，αzj，αzk）T，v为单元轴向变形增量，为单元扭转角增量，αyj，αyk分别为单元j端和k端绕单元坐标系ye轴的弯曲转角增量，αzj，αzk分别为单元j端和k端绕单元坐标系ze轴的弯曲转角增量；fV为单元基本内力增量向量，fV=（nx，tx，myj，myk，mzj，mzk）T，nx为单元轴向力增量，tx为单元扭矩增量，myj，myk分别为单元j端和k端绕单元坐标系ye轴的弯矩增量，mzj，mzk分别为单元j端和k端绕单元坐标系ze轴的弯矩增量；K为基本内力与基本变形关系的单元刚度矩阵；cy1，cy2，cz1，cz2均为梁单元弯曲刚度系数；Nx为单元轴向力（拉力为正）；l0为当前增量步单元初始长度；A为单元截面面积；Ix为单元扭转惯性矩；ρ为单元截面极回转半径；E，G分别为材料的弹性模量和剪切模量。

2.2 单元坐标系下单元杆端力向量与杆端位移向量的关系

将单元基本内力增量向量fV变换为单元坐标系下的单元杆端力增量向量，即

式中：fe为单元坐标系下单元杆端力增量向量，fe=（fjx，fjy，fjz，mjx，mjy，mjz，fkx，fky，fkz，mkx，mky，mkz）T，fjx，fjy，fyz分别为单元j端沿着单元坐标系xe，ye，ze轴方向的线力分量，mjx，mjy，mjz分别为单元j端绕单元坐标系xe，ye，ze轴的弯矩分量，fkx，fky，fkz，mkx，mky，mkz均为单元k端的相应分量；T为由单元基本内力向量变换为单元杆端力向量的变换矩阵。

2.3 单元坐标系下单元杆端位移与结构坐标下单元杆端位移的关系

设在上一增量步迭代结束时单元长度为l0，在当前增量步中单元长度为l1，则单元轴向变形增量为l1-l0，l1为单元j端和k端的整体坐标系下节点位移增量向量uj和uk的函数，即l1=l1（uj，uk）。设当前增量步中单元j端和k端的节点分别发生整体坐标系下的转角增量为θj和θk，将它们变换到单元坐标系下分别为rTθj和rTθk。应注意，θj和θk均为整体坐标系下的总量转角，它们包含着单元的弹性变形和刚体变位的贡献，其中单元的刚体变位在单元j端和k端产生的转角增量是相同的。设单元刚体变位在单元j端和k端产生的整体坐标系下的转角增量向量为β，将β变换为单元坐标系下为rTβ。从总转角增量rTθj和rTθk中减去刚体变位转角增量rTβ，剩余部分rT（θj-β）和rT（θk-β）则为单元弹性变形的转角增量。因此，当前增量步单元坐标系下单元弹性变形增量向量ue可写为

刚体变位转角增量向量β确定方法如下：设在上一增量步迭代结束时单元坐标系xe轴方向向量为r01，在当前增量步中单元坐标系xe轴方向向量为r1。由r01与r1所决定平面的法向量b为

按照先j端，后k端，先线位移，后角位移的顺序，将整体坐标系下杆端位移增量向量写成c=（uj，θj，uk，θk）T。根据矩阵微分理论，求fS关于杆端位移增量向量c的导数，则有

dfSdc=dSxdc[I12T（f0+KI67ue）]+SxTKI67duedc

（31）

式中：I12为12阶单位矩阵。

式（31）是一个关于c=（uj，θj，uk，θk）T的函数矩阵。令c=0，将其代入式（31）即可得到空间梁单元的切线刚度矩阵。

由于式（31）中dSxdc和duedc可实现精确求导，不需要忽略任何高阶项，因此，可得到空间梁单元切线刚度矩阵的精确计算公式。由于在式（31）求导过程时刚度矩阵K被视为常量而未涉及其求导，所以由式（31）确定的梁单元切线刚度矩阵也适用各种本构关系的梁单元，如弹性梁单元、塑性铰单元、精细塑性铰单元等。3 算例分析

文献[11]中采用有限元积分方法对六层空间刚架（图2）和二十层空间刚架（图3）进行弹塑性屈曲分析，得到结构荷载-位移曲线，其极限荷载因子为1.0。文献[10]，[12]，[13]中采用不同方法对该算例也进行了相应研究，并以文献[11]中的结果作为精确解进行对比。本文中为了验证所建立的梁单元切线刚度矩阵的有效性，对这2个算例也进行了弹塑性和几何非线性后屈曲分析。分析中荷载-位移平衡路径跟踪技术采用广义位移控制法（GDC）[14]，预测子为本文梁单元切线刚度矩阵，修正子为式（30），荷载因子初始步长λ0=0.01，收敛控制误差ε=0.001。弹塑性屈曲分析中，梁单元本构关系K[式（18）]采用精细塑性铰梁单元刚度矩阵[10，12]，几何非线性屈曲分析中，梁单元本构关系K采用具有稳定函数单元

刚度矩阵[15]。

3.1 六层空间刚架

六层空间刚架结构的尺寸及单元截面规格如图2所示。结构钢材屈服应力为250 MPa，弹性模量为206.85 GPa，剪切模量为79.293 GPa。楼面及屋面上均布重力荷载为9.6 kPa，按节点控制面积等效在每个楼层的柱顶。风荷载沿着y方向作用，等效在每个梁-柱节点的集中荷载为53.376 kN。按上述等效节点荷载分布比例进行加载，每个构件取为一个单元。endprint

对六层空间刚架进行弹塑性后屈曲分析，结构A点荷载-位移曲线如图4所示，极限荷载因子为1.005，与文献[10]～[13]中的结果非常接近，比精确解[11]仅高出0.5%。在极限荷载时A点x方向和y方向产生的位移分别为-0.129 78 m和0.223 86 m。分析结束时结构发生了很大的弹塑性扭转变形，且梁端和柱端形成了许多完全塑性铰，如图5所示，其中数字1～10为最先形成完全塑性铰的前10个位置。

3.2 二十层空间刚架

二十层空间刚架结构尺寸及单元截面规格如图3所示。结构钢材为A50钢，钢材的屈服应力fy=344.8 MPa，弹性模量E=200 GPa，剪切模量G=79 GPa。楼层及屋面均布重力荷载为4.8 kPa，按节点控制面积等效为集中荷载作用在每个楼层的柱顶上。结构立面上风荷载为0.96 kPa（沿y方向作用），按节点控制面积等效为集中荷载作用在梁-柱节点上。按等效节点荷载分布比例进行加载，每个构件取为一个单元。

对二十层空间刚架进行弹塑性后屈曲分析，A点和B点的荷载-位移曲线如图8所示。极限荷载因子为1.02，与文献[10]～[13]中的结果基本一致，比精确解[11]高出2%。在极限荷载处，A点和B点y方向所产生的位移分别为1.543 7 m和0.340 65 m。分析结束时结构发生了很大的弹塑性变形，并在梁端形成了许多完全塑性铰，如图9所示，其中数字1～10是最先形成完全塑性铰的前10个位置。

对二十层空间刚架进行几何非线性屈曲分析，结构A点和B点荷载-位移曲线如图10所示。极限荷载因子为2.409 4，在弹性极限荷载处，A点和B点y方向所产生的位移分别为3.026 3 m和1.039 5 m。分析结束时结构变形如图11所示。由图11可以看出，二十层空间刚架在下部1/3高度截面处，结构横向出现了明显鼓曲。

4 结 语

（1）本文中空间梁单元切线刚度矩阵分析时没有忽略任何高阶项，因此所得到的切线刚度矩阵具有足够精度，通过六层和二十层空间刚架的后屈曲分析得到验证。

（2）由于总体坐标系下杆端力与杆端位移关系式中的梁单元刚度矩阵K为常量矩阵，在应用矩阵微分理论求导时未涉及K的求导，因此本文中的空间梁单元切线刚度矩阵[式（31）]也适用其他的梁单元本构关系。

（3）从算例弹塑性后屈曲分析荷载-位移曲线可以看出，本文中的梁单元切线刚度矩阵可预测出结构极限荷载后较长距离的荷载-位移平衡路径。

（4）本文中的空间梁单元切线刚度矩阵可有效用于大型空间刚架结构的后屈曲分析。

参考文献：

References：

[1] GU J X，CHAN S L.New Tangent Stiffness Matrix for Geometrically Nonlinear Analysis of Space Frames[J].Journal of Southeast University：English Edition，2005，21（4）：480-485.

[2]陈常松，陈政清，颜东煌.势能增量驻值原理与切线刚度矩阵的解构规则[J].中南大学学报：自然科学版，2005，36（5）：892-898.

CHEN Chang-song，CHEN Zheng-qing，YAN Dong-huang.Principle of Stationary Incremental Potential Energy and Formation Rules of Tangent Stiffness Matrix[J].Journal of Central South University of Technology：Natural Science，2005，36（5）：892-898.

[3]陈 昕，王 娜.空间梁单元的切线刚度矩阵[J].哈尔滨建筑工程学院学报，1993，26（3）：24-29.

CHAN Xin，WANG Na.Tangent Stifness Matrix of Three Dimensional Beam Elements[J].Journal of Harbin Architecture and Civil Engineering Institute，1993，26（3）：24-29.

[4]YANG Y B，LIN S P，CHEN C S.Rigid Body Concept for Geometric Nonlinear Analysis of 3D Frames，Plates and Shells Based on the Updated Lagrangian Formulation[J].Computer Methods in Applied Mechanics and Engineering，2007，196（7）：1178-1192.

[5]陈军明，吴代华.三维梁单元的弹塑性切线刚度矩阵[J].华中理工大学学报，2000，28（2）：111-113.CHEN Jun-ming，WU Dai-hua.Elasto Plastic Tangent Stiffness Matrix of Three Dimensional Beam Elements[J].Journal of Huazhong University of Science and Technology，2000，28（2）：111-113.

[6]朱 文，周水兴.空间梁单元显式切线刚度矩阵推导[J].公路交通技术，2008（5）：40-49，53.

ZHU Wen，ZHOU Shui-xing.Derivation of Explicit Tangent Stiffness Matrix for Space Beam Element[J].Technology of Highway and Transport，2008（5）：40-49，53.endprint

[7]KIM S，RYU J，CHO M.Numerically Generated Tangent Stiffness Matrices Using the Complex Variable Derivative Method for Nonlinear Structural Analysis[J].Computer Methods in Applied Mechanics and Engineering，2011，200（1/2/3/4）：403-413.

[8]ORAN C.Tangent Stiffness in Space Frames[J].Jouranal of the Structural Division，1973，99（6）：973-985.

[9]KASSIMALI A，ABBASNIA R.Large Deformation Analysis of Elasticspace Frames[J]. Journal of Structural Engineering，1991，117（7）：2069-2087.

[10]刘永华.空间钢框架高等分析方法研究[D].哈尔滨：哈尔滨工业大学，2007.

LIU Yong-hua.Study on Advanced Analysis Method of Space Steel Frame[D].Harbin：Harbin Institute of Technology，2007.

[11]JIANG X M，CHEN H，LIEW J Y R.Spread-of-plasticity Analysis of Three-dimensional Steel Frames[J].Journal of Constructional Steel Research，2002，58（2）：193-212.

[12]刘永华，张耀春.空间钢框架精细塑性铰法高等分析[J].工程力学，2006，23（增1）：108-116.

LIU Yong-hua，ZHANG Yao-chun.Refined Plastic Hinge Analysis of Spatial Steel Frames[J].Engineering Mechanics，2006，23（S1）：108-116.

[13]CHIOREAN C G，BARSAN G M.Large Deflection Distributed Plasticity Analysis of 3D Steel Frameworks[J].Computers and Structures，2005，83（19/20）：1555-1571.

[14]YANG Y B，LIN S P，LEU L J.Solution Strategy and Rigid Element for Nonlinear Analysis of Elastically Structures Based on Updated Lagrangian Formulation[J].Engineering Structures，2007，29（6）：1189-1200.

[15]陈 骥.钢结构稳定：理论与设计[M].5版.北京：科学出版社，2011.

CHEN Ji.Stability of Steel Structures：Theory and Design[M].5th ed.Beijing：Science Press，2011.endprint