欧拉-伯努利梁单元刚度矩阵推导

张军锋，尹会娜，李 杰，陈 淮

(郑州大学 土木工程学院， 河南 郑州 450001)

梁理论有欧拉-伯努利(Euler-Bernoulli)梁和铁摩辛柯(Timoshenko)梁两种理论，后者计入剪切变形而前者一般不计入剪切变形[1]，前者又称为经典梁理论。现阶段常用的结构力学教程往往针对经典梁理论，并在矩阵位移法中通过两端固结等截面梁的梁端位移与梁端支反力的关系给出梁单元的刚度矩阵[2]。而在有限元理论中，则是通过对单元变形指定合理的形函数并结合最小势能原理[1，3]或虚功原理[4]给出单元刚度矩阵。尽管最小势能原理与虚功原理是等价的[3]，但推导方式仍有差异，尤其是最小势能原理涉及变分原理而不便理解。另外，部分文献并未给出完整坐标系统以及荷载、内力和位移的正方向，不同文献所给坐标系统也不完全一致，这也不便理解其推导过程。

为系统介绍等截面梁刚度矩阵的推导过程，基于欧拉-伯努利梁理论，采用统一的坐标系统和参数正方向规定，详细介绍了不同受力模式包括伸缩、扭转和弯曲问题的拉格朗日形函数或厄米特形函数，系统给出了刚度矩阵推导所需的基本方程包括几何方程、本构方程和平衡方程，并最终得到单元刚度矩阵，以便深刻全面的理解梁单元刚度矩阵以及对梁结构受力特性的认识，并有助于参阅其他专业书籍。

1 坐标系统与形函数

欧拉-伯努利梁理论基于Kirchhoff假设，即假定变形前垂直于中面的截面，变形后仍为平面且与中面保持垂直，或者说可忽略剪切变形的影响，适用于高度远小于跨度的梁[1]。

以xy平面梁单元为例，图1给出了右手螺旋法则的坐标系统，并针对伸缩、扭转和弯曲问题给出了对应的杆端荷载和分布荷载，所示荷载方向均与坐标轴正方向一致，也即荷载的正方向。图1(c)还给出了截面内力正方向：轴力F以受拉为正，扭矩T的正方向与F方向一致，弯矩M以梁底(即-y方向)纤维受拉为正，剪力Q以使微段发生逆时针转动为正。部分文献对单元刚度矩阵推导时，单元长度范围内还受有集中荷载，包括集中力和集中力矩，而实际对结构划分单元时往往在集中荷载位置均设置有节点，故图1中的单元没有再受集中力作用。另外，单元范围内的荷载也只影响最终的荷载向量而对刚度矩阵没有影响[1，3-4]。

图1单元坐标系以及参量正方向示意图

形函数是用来表示结构变形形状的物理量，通过结点位移与结点之间单元的变形来表示整个结构的变形。形函数实际上就是数学上的插值函数，只是在有限元中又称形函数。常用的形函数有拉格朗日(Lagrange)形函数和厄米特(Hermite)形函数，众多有限元书籍以及ANSYS软件中的Beam 4单元对伸缩和扭转问题采用拉格朗日形函数，而对弯曲问题采用厄米特形函数。这也是因为这三种受力模式是非耦合的，故可以独立进行分析。对不同受力问题差值函数的构造过程可见文献[3]。

对于所示x-y平面内的两节点单元，其位移函数φ采用拉格朗日形函数表达式为：

(1)

(2)

(3)

(4)

(5)

Ni(xj)=Ni(sj)=δij(i,j=1,2)；

(6)

其中：i表示结点编号1和2(或I和J)；φi表示i结点的位移，可表示轴向位移u或扭转角θ；s为单元相对位置并取-1≤s≤1(见图1、式(3))；L为单元长度；xc为单元中点坐标。Ni(s)和Ni(x)表示拉格朗日形函数，且应满足式(6)，xj和sj为单元j节点的绝对坐标和相对坐标，其中δij为Kronecker函数，并且式(6)中的前者是拉格朗日形函数本身的性质，后者则是差值函数完备性的要求[1]。

另外，Ni(s)和Ni(x)本质是相同的，只是表现形式不同(见式(4)、式(5))：如将式(5)中的相对坐标s替换为式(3)的关于x的函数，则可得式(4)即为关于杆件绝对坐标x的形函数，下文的厄米特形函数也同样如此。下文除对轴向位移给出两种形式的位移函数外(见式(7))，其他均以相对坐标s给出(见式(5)、图2)。另外，本文采用的s取值方法也是ANSYS技术手册中采用的方式，但在部分文献[1]中也常采用0≤s=x/L≤1的形式，此时对应的形函数表达式会有变化，但并不影响最终所得单元刚度矩阵，在参阅时需注意。

图2形函数示意图

式(4)、式(5)所示拉格朗日差值函数适用于单元伸缩和扭转问题，由式(1)—式(5)可得轴向位移u和扭转角θx的位移函数分别如式(7)和式(8)所示：

(7)

(8)

两节点单元弯曲问题的竖向位移函数φ则采用厄米特形函数表示为：

(9)

(10)

(11)

(12)

(13)

为与节点自由度表示形式一致，对式(7)通常表示为：

(14)

(15)

实际上，式(7)、式(8)和式(15)也正是ANSYS技术手册中所给的Beam 4单元位移函数表达式。为便于刚度矩阵分析，上述三式可统一用向量形式表达为：

(16)

其中：N和Φ分别为形函数向量和节点位移列向量；ψr表示节点I和J的位移，可分别为轴向位移u或扭转角θ、竖向位移ω和转角α；n为独立节点位移数量；对于三类问题的各自取值分别为式(17)、式(18)、式(19)，同时还给出了形函数向量N相应的一阶和二阶导数N′、N″，并且其各阶导数依然为向量。

(17)

(18)

(19)

2 刚度矩阵

由于梁单元的三种受力模式是非耦合的，故可以分别求得单一受力模式下的单元刚度矩阵，再综合为任意荷载模式下的单元刚度矩阵。不同文献对基于形函数的单元刚度矩阵推导过程并不一致，比如文献[1]利用最小势能原理进行推导，文献[3-4]还利用虚功原理进行了推导，但最终所得刚度矩阵的计算方法和具体表达式(表1中刚度矩阵计算方法栏和刚度矩阵栏)是一致的。文献[3-4]还介绍了两种推导方法的等价性：两者本质上是一致的，但前者可以判断解的收敛性以及解的下限性质。由于推导过程过于复杂且在文献[1,3-4]有较为详细的过程，故下文不再给出，仅以虚功原理为例，基于图1所示坐标系和参数正方向，补充给出三种受力模式的几何方程、本构方程、截面内力、平衡方程、刚度矩阵计算过程及结果。

表1 Beam 4单元三种受力模式下的基本方程

注：曲率κ本身应为正值，但为适应应变计算还是对其给出了符号。整个表格中的N和B均为绝对坐标x的函数，N′、N″为相对坐标s的表达式。刚度矩阵计算过程一栏中，倒数第二个等式仅针对等截面杆件；最后一个等式对不同相对坐标s的形式(即-1≤s≤1或0≤s≤1)所得表达式会有差异，但不影响最终所得刚度矩阵。

表1中：u、θ、ω分别表示轴向位移、扭转角和挠曲位移；β和α分别为扭矩作用下单位长度的扭转角和弯矩作用下的截面转角；ε和σ分别为正应变和正应力，下标x和b分别表示轴力和弯矩作用；γ和τ分别为剪应变和剪应力，下标t和b分别表示扭矩和弯矩作用；E和G为弹性模量和剪切模量，且有式(20)，μ为泊松比；A、J、I表示截面面积、极惯性矩和抗弯惯性矩；y为偏离截面中性轴的坐标值；ρ为距离截面扭心的距离。Fe和Ke分别为对应受力模式的荷载向量和单元刚度矩阵，节点位移列向量Φ同样是在单元坐标系下的位移。N、B分别为形函数向量和应变转换向量，并且对不同受力模式的取值是不同的(见式(17)、式(18)、式(19))，只是为使形式统一而采用了相同的符号。对于各向同性材料的一维线单元，应力与应变之间的弹性矩阵D对正应力和剪应力分别为E和G；如果是二维壳单元或三维实体单元，则D将是一个矩阵[1]。

(20)

另外，最终根据形函数所得刚度矩阵与结构力学方法所得刚度矩阵完全一致，这也说明用形函数(式(17)、式(18)、式(19))表示的位移函数φ是真实的，满足理论解。当然，如果所用形函数是近似的，所得刚度矩阵也与结构力学方法所得刚度矩阵不一致：对于弯曲问题如果同样采用拉格朗日形函数就会出现这样的结果。鉴于三种受力模式的独立性，将表1所列各自刚度矩阵综合在一起则可得式(21)，其对应的总体刚度矩阵、节点荷载向量和位移向量如式(22)所示。

Fe=KeΦe

(21)

Fe=[FxI,FyI,MI,TI,FxJ,FyJ,MJ,TJ]T

(22)

Φe=[uI,ωI,αI,θI,uJ,ωJ,αJ,θJ]T

3 结 论

梁单元刚度矩阵是结构力学中的一个重要内容，本文对等截面梁，基于欧拉-伯努利梁理论，采用统一的坐标系统和参数正方向规定，详细介绍了各种受力模式的形函数，系统给出了等截面梁的刚度矩阵推导过程，为配合有限元专业文献的理解提供了参考。