变截面梁单元刚度矩阵的推导及影响因素分析

陈代海， 周帅， 李银鑫， 刘明杰， 李波

(1 郑州大学 土木工程学院，河南，郑州 450001; 2.中铁第四勘察设计院集团有限公司，湖北 武汉 430063)

变截面结构由于其能合理分配自重，充分发挥材料强度优势，受力合理，因此在实际工程中应用广泛。在对其进行有限元数值分析时，需要推导变截面梁单元刚度矩阵，王勖成[1]、叶康生等[2]进行了理论基础研究；传光红等[3]等采用势能驻值的原理，推导了变截面Timoshenko梁的单元刚度矩阵；Al-Gahtani H J[4]研究了变截面Bernoulli-Euler梁的轴力、弯曲及扭转单元刚度矩阵；张琪等[5]通过Hamilton原理，推导出变截面梁段单元的自由振动方程，得到了单元刚度矩阵；李爽等[6]利用平衡条件得到梁截面内力的分布模式，在其变形前建立平衡方程，再由虚力原理导得刚度矩阵；马志贵等[7]假定变截面的形式为梁高线性变化，根据有限元基本原理，推导出变截面梁单元的单元刚度矩阵公式；John S,et al[8]讨论了一种新型的变截面单元的刚度矩阵，并对其进行试验研究；陆念力[9]、武芳[10]中采用直接法或传递矩阵法推导单元刚度矩阵；传光红、陆念力、耿文宾通过将变截面梁单元进行简化，如单元两端截面等效为单元中间截面，将截面两端惯性矩做平均化处理等，从而得出变截面梁单元的刚度矩阵[3,9,11]；Al-Gahtani H J、李爽、John S、武芳没有给出变截面梁单元刚度矩阵的具体表达式或需要确定的参数较多[4,6,8,10]，不便于实际应用。因此，该文基于能量变分原理，推导变截面空间梁单元刚度矩阵，得到2结点与3结点变截面梁单元刚度矩阵的具体表达式，并通过对比相应Ansys单元类型的计算结果验证其正确性。在此基础上，采用Ansys建立变截面连续梁模型，分析变截面单元简化为等截面单元、梁高变化线形及变截面单元划分长度对变截面连续梁静力位移的影响。

1 变截面梁单元刚度矩阵的推导

有限单元法中，常采用2结点单元或3结点单元离散变截面梁。该文基于能量变分原理，建立轴力单元、扭转单元及弯曲单元的泛函，根据弹性力学基本方程及边界条件，求泛函的驻值，得到变截面梁单元的刚度矩阵表达式。

1.1 2结点变截面梁单元

对变截面梁单元刚度矩阵的推导，可在局部坐标中考虑梁单元承受的横向荷载及弯矩。以承受均布荷载q(x)和集中荷载Pj作用的悬臂梁为例，如图1所示，其挠度函数为ω(x)。

根据弹性力学基本方程可得：

(1)

式中:k为梁中面变形后的曲率；M为截面承受的弯矩；Q为截面承受的剪力；θ为端部给定的转角；It为截面惯性矩；E为弹性模量；x为梁端到计算截面的距离。

图1 弯曲变形梁单元示意图

根据能量变分原理建立泛函如下：

(2)

由有限单元法计算梁单元的弯曲刚度矩阵时，采用一维Hermite插值多项式表示挠度插值函数为：

(3)

式中：N1=1-3ξ2+2ξ3，N2=(ξ-2ξ2+ξ3)l，N3=3ξ2-2ξ3，N4=(ξ3-ξ2)l；l为悬臂梁长度；αi为结点位移；N为形函数；αe为αi的矩阵表示形式。

将式(3)代入式(2)，求解δ∏p=0，得到Kα=P，即可求得任意变截面梁单元的弯曲刚度矩阵Kω为：

(4)

推导变截面梁单元的轴力刚度矩阵及扭转刚度矩阵时，采用拉格朗日插值多项式表示其位移函数为：

(5)

采用类似于梁单元弯曲刚度矩阵的推导方法，可得到变截面梁单元的轴力刚度矩阵Kα及扭转刚度矩阵Kb。

以矩形变截面梁单元为例，其宽度为b，假设截面高度h=h0+mx，其中h0为起始端单元高度，m为梁高线性变化的比例系数。则沿单元局部坐标轴y轴和z轴的惯性矩分别为Iy=hb3/12和Iz=bh3/12，单元截面面积为A=bh。其轴力刚度矩阵Kα及扭转刚度矩阵Kb为：

(6)

(7)

式中:β为截面扭转系数;G为剪切模量;则单元刚度矩阵Ke可以表示为：

(8)

在实际工程中，受到弯曲变形的影响，梁单元所受横向剪切力将引起其附加挠度，即变形后的梁单元会发生翘曲。式(9)中并未考虑剪切变形的影响，且构造位移函数时，截面转动角度为：θ=dω/dx，若在此基础上考虑剪切变形，其推导过程将较为复杂。因此，在保证精度的前提下简化推导过程，对挠度ω与截面转动

(9)

角度θ进行独立插值，推导变截面梁考虑剪切变形时的单元刚度矩阵。

(10)

对式(10)求驻值，可得到梁单元弯剪刚度矩阵为：Kω=Kωy(z)+Ks。

以矩形变截面梁单元为例，参数同上，推导2结点变截面梁单元的刚度矩阵。其弯曲刚度矩阵Kωy(z)及剪切刚度矩阵Ks表示如下：

(11)

(12)

(13)

式中:k为横截面剪应力影响系数，对于矩形截面，k=6/5。则考虑剪切变形的2结点变截面梁单元刚度矩阵如式(9)所示。

1.2 3结点变截面梁单元

3结点梁单元刚度矩阵的推导方法与2结点单元类似，其位移函数采用拉格朗日插值多项式可表示为：

(14)

3结点矩形变截面单元的轴力刚度矩阵及扭转单元刚度矩阵可表示为：

(15)

(16)

单元弯剪刚度矩阵为：

(17)

(18)

(19)

则3结点变截面梁单元刚度矩阵可以表示为[见式(20)]：

2 变截面梁单元刚度矩阵的验证

为验证上述2结点单元与3结点单元刚度矩阵表达式的正确性，以如图2所示的变截面悬臂梁为研究对象，运用Ansys有限元软件，分别采用Beam44和Beam189单元建立其有限元模型。悬臂梁长l为1 m，i端截面高hi为0.2 m，宽bi为0.1 m；j端截面高hj为0.3 m，宽bj为0.1 m，弹性模量E为2×108Pa，泊松比μ为0.2，剪切模量G为8.333×107Pa。悬臂梁一端固定，一端自由，自由端施加荷载Fx=50 N，Fy=50 N，Fz=50 N，Mx=50 N·m，My=50 N·m，Mz=50 N·m。

基于上述推导过程得到单元刚度矩阵K[计算式见式(9)、(20)]，结合荷载列阵P，求解静力平衡方程Kα=P，得到悬臂梁自由端的位移，并与Ansys计算结果进行对比，如表1所示。

由表1可以看出：根据能量变分法推导的变截面梁单元刚度矩阵，对于2结点单元与3结点单元，该文计算结果与Ansys计算结果基本一致，验证了该文变截面梁单元刚度矩阵表达式的正确性，其可用于变截面结构的有限元分析。

3 变截面梁单元在数值模拟中的影响因素分析

在连续梁桥的有限元数值模拟中，变截面梁单元应用较多。以往的有限元分析中，为简化运算，常将变截面梁单元换算为等截面单元，由此带来的计算误差值得研究。另外，由上述变截面梁单元刚度矩阵的表达式可以看出，单元长度和截面高度是影响梁单元刚度的主要因素。因此，该文以3跨变截面连续梁桥为研究对象，采用Ansys软件建立其有限元模型，分析变截面单元简化为等截面单元、梁高变化线形及变截面单元划分长度等因素对桥梁结构静位移的影响。

(20)

图2 变截面悬臂梁示意图(单位：m)

表1 悬臂梁自由端的位移对比

3.1 工程实例

某变截面连续梁桥跨度为：(40+64+40) m，桥梁宽7.1 m，变截面段长29.5 m，桥梁纵向布置如图3所示。边支点处截面高2.8 m，中支点处截面高5.2 m，截面尺寸如图4所示。采用C50混凝土，其弹性模量E为3.45×1010Pa，泊松比μ为0.2，材料重度γ为25 kN/m3。考虑桥梁结构自重和均布荷载Q=100 kN/m的作用。

图3 变截面连续梁纵向示意图(单位：cm)

图4 变截面连续梁截面示意图(单位：cm)

3.2 变截面单元简化为等截面单元的影响

基于上述工程实例，变截面梁高按2次抛物线变化，对于每个变截面单元，单元长度分别选取为3、4、5和6 m，采用两种等截面分别对变截面进行简化：一种等截面定义为等截面1，为变截面单元的中部截面；另一种等截面定义为等截面2，其截面特性在有限元模型中以实常数的形式表示，实常数中的截面面积取两端截面面积的算术平均值，惯性矩取两端截面的几何平均值。用Ansys软件中的2结点变截面单元Beam44离散变截面连续梁，建立变截面单元简化为等截面单元的连续梁有限元模型。在均布荷载作用下，计算桥梁的竖向位移，结果如表2所示。

表2 不同计算模型的桥梁竖向位移对比

由表2可以看出：两种等截面模型的计算结果相近，说明两种等截面均可用于简化变截面。变截面单元简化为等截面单元时，当变截面单元划分长度为3、4 m，即为变截面段长度的1/10和1/7时，等截面与变截面模型的位移值相差较小，最大相对差值小于3%；当变截面单元划分长度为5、6 m，即为变截面段长度的1/6和1/5时，等截面与变截面模型的最大位移差值较大，最大相对差值为13.48%。因此，在建立变截面梁单元模型时，当单元划分长度达到变截面段长度的1/6～1/5时，不宜采用等截面对其进行简化。

3.3 梁高变化线形的影响

变截面连续梁梁高变化线形常以抛物线为主，为研究抛物线形状引起的梁高变化对桥梁刚度的影响，采用Ansys软件中的3结点单元Beam189离散变截面连续梁，建立梁高变化线形分别为1.5、1.8和2次抛物线时的3种有限元模型，在均布荷载作用下，计算桥梁的竖向位移，在自重作用下，计算桥梁支座反力。不同梁高变化线形下的桥梁竖向位移对比如图5所示，支座反力对比如表3所示。

图5 不同梁高变化线形下的桥梁竖向位移对比

表3 不同梁高变化线形下的桥梁支座反力对比

由图5可知：在均布荷载作用下，桥梁跨中最大位移值随着梁高变化抛物线形阶次的增加而增大，最大相对差值为6.45%。由表3可知：在自重作用下，桥梁支座反力随着梁高变化抛物线形阶次的增加而减小，最大相对差值为2.5%。说明当梁高变化线形分别为1.5、1.8和2次抛物线时，梁高变化线形对变截面连续梁的质量分布影响较小，而对其刚度分布影响较大。

3.4 变截面单元长度的影响

为分析变截面单元划分长度对桥梁刚度的影响，采用Ansys中的3结点单元Beam189离散变截面连续梁，建立单元长度分别为3、4、5和6 m时的4种有限元模型，计算均布荷载作用下的桥梁竖向位移，得到其竖向位移对比如图6所示。

图6 不同变截面单元划分长度下的桥梁竖向位移

从图6可以看出：当变截面单元划分长度为3、4 m，即为变截面段长度的1/10和1/7时，两种模型的位移值基本吻合；当变截面单元划分长度为5、6 m，即为变截面段长度的1/6和1/5时，桥梁跨中的最大竖向位移差值较大。当变截面单元划分长度大于变截面段长度的1/7后，随着单元长度的增大，桥梁跨中竖向位移最大值逐步减小。因此，为保证计算精度和效率，在建立变截面梁有限元模型时，宜将变截面单元划分长度控制在变截面段长度的1/7以内。

4 结论

(1) 基于能量变分法，推导了2结点和3结点变截面梁单元的刚度矩阵具体表达式。并通过实例计算，对比分析Ansys软件计算结果，验证了其正确性，说明其可用于变截面结构的有限元分析，可为有限元计算程序的编制提供参考。

(2) 在建立变截面连续梁有限元模型时，宜将变截面单元划分长度控制在变截面段长度的1/7以内。当单元划分长度达到变截面段长度的1/6～1/5时，不宜采用等截面对其进行简化。

(3) 当梁高变化线形分别为1.5、1.8和2次抛物线时，梁高变化线形对变截面连续梁的质量分布影响较小，而对其刚度分布影响较大。