考虑剪切变形的变截面欧拉梁单元刚度矩阵

张军锋 李 杰 尹会娜 陈 淮 叶雨山

（1.郑州大学土木工程学院，郑州450001；2.中国建筑第七工程局有限公司，郑州450001）

0 引 言

欧拉梁理论又称为欧拉-伯努利（Euler-Bernoulli）梁理论或经典梁理论，其忽略剪切变形对挠度的影响，但对高度较大的梁，其剪切变形对挠度的贡献将不可忽略［1-2］。为提高欧拉梁理论的适用性，有限元中往往基于欧拉梁理论并计入剪切变形的影响，比如ANSYS中的Beam4∕44单元，前者针对等截面梁单元而后者针对变截面梁单元。

等截面梁单元刚度矩阵的推导可采用结构力学形常数方法［2］和以形函数为基础的最小势能原理或虚功原理方法［1，3-4］，但对变截面梁单元刚度矩阵的推导却鲜有涉及。文献［5］针对两结点梁单元给出了变截面梁单元的单刚矩阵，但未明确刚度矩阵的具体推导原理和过程，也没有计入剪切变形，且所给单刚矩阵仅针对弯曲问题而不涉及伸缩和扭转问题。文献［6］采用传递矩阵法推导了任意变截面梁单元的刚度矩阵，但此方法不涉及形函数，与通用有限元理论不一致。文献［7］则以形函数为基础，借助于MATLAB的符号运算功能建立了梁高成线性或抛物线形变化的变截面梁单元刚度矩阵，但其矩阵元素颇为复杂，不便使用。ANSYS的Beam44与Beam4的单刚矩阵形式相同，只是前者对截面参数包括面积A和惯性矩I采用了“平均值”，但并未给出截面参数“平均值”的推导原理。

基于有限元理论中的形函数和最小势能原理，在等截面欧拉梁单元刚度矩阵推导过程的基础上［8］，系统阐述了计入剪切变形的变截面欧拉梁单刚矩阵推导过程，并经过对比分析明确了ANSYS中Beam44单元刚度矩阵的推导方法。为便于阐述，仅针对平面梁单元进行推导。另外，由于ANSYS中的Beam44这一变截面梁单元要求截面在两节点之间是线性梯度变化的，下文建立的变截面梁同样符合这一要求。

1 基本方程

以XY平面梁单元为例，图1给出了右手螺旋法则的坐标系统，并针对伸缩、扭转和弯曲问题给出了对应的杆端荷载和分布荷载，所示荷载方向均与坐标轴正方向一致，也即荷载的正方向。图1（c）还给出了截面内力正方向：轴力F以受拉为正，扭矩T的正方向与F方向一致，弯矩M以梁底（即-y方向）纤维受拉为正，剪力Q以使微段发生逆时针转动为正［1，3-4］。

图1 单元坐标系以及参量正方向示意图Fig.1 Positive direction regulations in element coordinate system

对于所示X-Y平面内的两节点单元，其位移函数φ也即变形包括伸缩、扭转、弯曲和剪切位移四种变形且均可采用形函数表达为

式中：N和Φ分别为形函数向量和节点位移列向量；ψr表示节点i和j的位移，可分别为节点的轴向位移u或扭转角θ、弯曲竖向位移ωb和剪切竖向位移ωs和截面转角α；n为独立节点位移数量，除弯曲竖向变形ωb的r=4外，其他三种变形均为r=2（表 1）；s为单元相对位置并取-1≤s≤1（图 1，式（2）），xj和sj为单元j节点的绝对坐标和相对坐标；L为单元长度；xc为单元中点坐标；Ni（s）和Ni（x）表示拉格朗日形函数或厄米特形函数，表1给出两种形函数的形式和适用变形［1］。

需要指出，从结构力学分析可知，变截面梁的形函数极为复杂，故在有限元中对变截面梁依然直接使用等截面梁的形函数以简化分析。

表1 不同变形的形函数Table 1 Shape functions for different deformation functions

需要说明，剪切变形的计入使梁单元看似存在4种受力模式，但其纯剪切和纯弯曲受力实际上是耦合的［1］。为便于阐述，下文依然称为4种受力模式，并对每种受力模式分别分析其单刚矩阵，再综合为任意荷载模式下的单刚矩阵。为节约篇幅，基于前述位移函数和形函数以及最小势能原理，以图1所示坐标系和参数正方向，表2直接给出各种受力模式的单刚矩阵计算方法，最小势能原理以及所涉及的各种受力模式的几何方程、本构方程、截面内力和平衡方程详见文献［8］；纯剪切变形刚度矩阵的推导中，还根据能量等效原则引入了剪应力不均匀修正系数k（以下简称剪切系数k）以便使用统一的单刚矩阵计算公式，详见文献［1］。

表2中，E和G为弹性模量和剪切模量，μ为泊松比；截面参数A、J、I分别为面积、抗扭惯性矩和抗弯惯性矩；由于仅针对平面梁单元，此处的I为绕z轴的抗弯惯性矩也即Iz（图1、图2）；y为偏离截面中性轴的坐标值；S（y）和b（y）为截面剪应力计算位置一侧截面积对中性轴的静矩和计算位置的截面宽度。对于变截面梁单元，截面参数A、J、I均随s变化，故在刚度矩阵推导过程中应始终将其放在积分号内；对于截面剪切系数k，从其计算公式可知，其仅与截面形状有关，如矩形、圆形和薄壁圆环截面分别为6∕5、10∕9和2，但部分变截面梁的k将随s位置变化（图2），故亦置于积分号内。当然，如果为等截面梁，截面参数A、J、I和剪切系数k均可置于积分号外，从而得对应的刚度矩阵如表2所列。下文则以矩形、圆形和薄壁圆环和薄壁箱形截面为例，分析变截面梁的刚度矩阵，并与ANSYS所给刚度矩阵进行对比。

2 刚度矩阵中的截面参数分析

2.1 ANSYS中Beam44单元截面参数取值

在ANSYS中，基于欧拉梁理论提供的等∕变截面梁单元分别为Beam4和Beam44，两者的单刚矩阵形式相同（式（3）、式（4）），只是Beam44刚度矩阵中的截面参数A、I、J均采用了“平均值”作为等效截面参数（式（5）），其中的参数下标i和j分别表示左右节点。或者说ANSYS的Beam44本质上也是一种等截面梁，只是以左右两端截面参数的“平均值”作为截面参数取值。在ANSYS中，Beam44单元只能通过实常数命令直接定义i和j节点的截面参数A、I、J，并对k只能取一个值，而无法通过命令SECTYPE和SECDATA定义i和j节点的截面形状和尺寸。这也说明，Beam44单元的单刚矩阵仅由两端截面参数和式（3）-式（5）决定，与截面形状沿杆件的变化规律无关，但从下文的理论分析可知，这样所得的单刚矩阵对大部分截面只是一种近似。

表2 不同受力模式的单刚矩阵推导Table 2 Deviation of element stiffness matrixes

2.2 剪切系数k

对于截面剪切系数k，ANSYS中的Beam44单元和Beam4单元均只能输入一个参数，即对变截面梁认为其不随位置改变。实际上，对于变截面梁，矩形、圆形和薄壁圆环截面的k本身并不随位置变化，但箱形截面、工形和T形截面的k随位置变化：图2依文献［10］对一箱形截面梁给出了不同截面尺寸的k值。但是，在单元两端截面参数A、J和I比值0.5～2.0的约束下（详见下文），k的变化较为有限，故下文亦对k取常数即取单元中截面的k值，从而可将表2中纯剪切刚度矩阵推导中的k置于积分号以外，同时对Beam44单元亦有式（4）。

图2 箱型截面的截面参数Fig.2 Section parameters of a box section

2.3 矩形截面参数分析

对截面剪切系数k取定值后，表2中各受力模式刚度矩阵的积分式内仅有截面参数和形函数导数向量，但因截面形状多样，且截面尺寸随杆件长度变化形式多样，直接对其进行积分仍非常复杂，且结果与ANSYS中Beam44单元所得单刚矩阵中的部分元素仍有较大差异。下文首先以最简单的矩形截面为例进行分析，并且此时剪切系数k本身不随s变化。

取矩形截面变截面梁i、j两端截面尺寸分别为bj×hj=c1bi×c2hi（式（6）），其中c1和c2分别为截面宽度和高度的变化系数（以下简称截面变化系数），假定截面尺寸随s线性变化（这也与ANSYS中对Beam44单元的假定一致），从而可得各截面参数如式（7）所示。

对于式（7）需要说明：对于本文的平面梁单元，其抗弯惯性矩I也即Iz，但同时也给出了Iy以便对比；矩形截面抗扭惯性矩J并没有通用的计算表达式，往往采用J=χb4计算，其中系数χ表达式极为复杂，故实际计算时往往根据高宽比h∕b直接对χ取值［11］。由于难以给出χ随s的变化表达式，故式（7）对J的表达式仅适用于截面变化系数c1=c2，也即截面高宽比h∕b不随s变化的情况。

将式（7）所示A、J、I带入表2进行积分计算即可得单一受力模式下的单刚矩阵。从表2所列刚度矩阵推导过程可知，对于除纯弯曲以外的其他3种受力模式，积分计算仅针对截面面积A和抗扭惯矩J。对于轴向伸缩和纯剪切受力，积分计算仅针对截面面积A且可得式（8）。对比式（5）可知，只有当截面变化系数c1=c2时也即截面宽度和高度按同一规律变化时，ANSYS所给等效截面面积才与理论值一致，在c1≠c2时ANSYS所给等效截面面积只是理论值的一种近似。对于扭转受力，由于难以给出c1≠c2时J的计算表达式，但c1=c2时则有式（9）成立，同样说明只有当c1=c2时ANSYS所给等效截面抗扭惯性矩（式（5））才与理论值一致。

纯弯受力的积分表达式更为复杂，其抗弯惯性矩I和形函数二阶导数N''均有积分变量s，积分所得表达式过于复杂，且通过计算对比发现即使c1=c2时所得刚度矩阵理论值也与ANSYS所得不一致。为明确ANSYS中Beam44单元单刚矩阵的计算方法，通过尝试发现，如果将分进行计算（式（10）），也即对几何矩阵和截面参数分别积分，则在c1=c2时所得刚度矩阵将与ANSYS结果一致（式（11）-式与ANSYS所给等效截面抗弯惯性矩（式（5））一致且所得矩阵与纯弯受力模式的刚度矩阵形式一致；两者联合并考虑系数8E∕L3后即得ANSYS变截面梁纯弯刚度矩阵。这实际也明确了ANSYS中对Beam44单元纯弯模式刚度矩阵的计算方法，即将″ds进行计算，但对矩形截面也仅当c1=c2时所得单刚矩阵与ANSYS所给单刚矩阵一致。

明确了矩形截面单一受力模式下ANSYS中单刚矩阵后，即可将其综合为式（3）所示任意受力模式下的单刚矩阵。而从上述推导可知，对于矩形截面梁，各受力模式下ANSYS中Beam44单元所得单刚矩阵仅在截面宽度和高度按同一规律变化也即截面变化系数c1=c2时与理论值一致，在c1≠c2时仅为一种近似。实际上，ANSYS中Beam44单元对矩形截面并无c1=c2的限制，但ANSYS要求Beam44 单元两端截面参数的比值 Ai∕Aj、Ji∕Jj、Ii∕Ij均应在0.5～2.0之间，否则会提示警告信息，如果上述比值超出0.1～10.0的范围，则会提示错误信息。这可能也是为了保证当c1≠c2时ANSYS所得等效截面参数与理论结果的差异不至过大而提出的要求：两端截面参数比值越接近，两者的差异也就越小。

2.4 其他常用截面参数分析

除矩形截面外，圆形、圆环和箱形截面也是常用的截面形式。对于圆形截面，因其截面尺寸仅有1个参数即半径R，截面变化系数同样仅有1个即c（式（14））：因圆形截面的Iz=Iy=0.5J=I故仅给出I的表达式，并且对比式（7）可知，只需将式（7）的c1和c2取为c即可的圆形截面A、J、I的表达式，且圆截面的k值亦不随半径变化。对其分析结果表明，在采用式（10）后，其理论推导所得刚度矩阵与ANSYS采用等效截面参数的表达式完全一致。

对于薄壁圆环截面，其截面尺寸有2个参数即外径R与壁厚t，相应的截面变化系数同样为2个，即 c1和 c2。故其截面尺寸和 A、J、I的表达式（式（15））与式（7）一致，并且由于是圆环截面，其J值表达式不再有c1=c2的限制。对于箱形截面，其截面尺寸参数更多，理论上有6个，但一般顶底板厚度和左右腹板厚度分别相等故有4个参数（图2），相应的截面变化系数同样为4个。对圆环和箱形其分析结果表明，在假定k值不随杆件位置变化且使用式（10）后，其理论推导所得刚度矩阵同样在c1=c2=c和c1=c2=c3=c4=c时与ANSYS采用等效截面参数的表达式一致。

需要说明，式（15）对A、J、I的表达式仅适用于薄壁圆环，并非精确表达式（式（16））。如果采用式（16）的精确表达式，则在c1≠c2时对A和I积分不再有式（8）和式（11）的结果，而式（15）的近似表达式在c1≠c2时同样有式（8）和式（11）成立，但精确表达式在c1=c2=c时同样可得式（5）。

综上可推知，ANSYS对Beam44变截面梁单元单刚矩阵推导时：对纯弯模式采用进行计算，同时对截面剪切系数k取定值；所给等效截面参数A、J和I表达式源于截面各个尺寸的截面变化系数相同时的情况，或者梁单元两端截面形状为相似形且截面尺寸随杆件线性变化的情况，其他情况ANSYS所得等效截面参数与理论值会有一定偏差。

3 结 论

基于有限元理论系统阐述了计入剪切变形的变截面欧拉梁单刚矩阵推导过程，并以矩形、圆形、圆环和箱形截面梁为例与ANSYS中Beam44单元刚度矩阵进行了对比，明确了Beam44单元刚度矩阵的推导方法。研究发现：

（1）变截面梁因截面形状多样且截面尺寸随杆件长度变化形式多样，即使对于线性梯度变化的变截面梁，对部分截面也难以给出截面参数，如抗扭刚度J和剪切系数k随杆件位置的表达式，且变截面梁使用等截面梁的形函数本身也是一种近似，故难以对其进行积分计算获得普适性的刚度矩阵理论表达式。

（2）ANSYS对Beam44变截面梁单元单刚矩阵推导时，纯弯模式的积分计算中对几何矩阵和截面参数分别积分，即以代替进行计算，这样既可简化计算，又可使各种受力模式的计算方式一致。

（3）Beam44单元对截面剪切系数k取定值以简化计算，故在使用中宜取单元中截面的剪切系数k或左右两端截面的均值以减少误差。

（4）ANSYS所得Beam44单元的单刚矩阵形式与等截面梁完全一致，只是对截面面积A、抗扭惯矩J和抗弯惯矩I采用了等效值。ANSYS所给等效截面参数A、J和I表达式，源于梁单元两端截面形状为相似形且截面尺寸随杆件线性梯度变化的情况，其他情况的等效截面参数与理论值会有一定偏差。ANSYS对Beam44单元仅要求截面在两节点之间按线性梯度变化，并无要求两端截面形状为相似形，这实际也是将简单模式所得表达式直接推广应用至各种情况。但ANSYS对Beam44单元要求两端截面参数的比率尽可能接近1.0，这可能也是为了控制等效截面参数与理论值的偏差，故对变截面梁在建模时应尽量减小单元长度以提高计算精度。