摘 要:固定形状的单元位移插值函数不能合理地近似变截面梁内部的位移变化,从而影响了传统梁单元用于计算变截面梁的精度。采用直接基于单元平衡的思想给出了计算变截面梁反应的有限元方法,解决了单元位移插值函数局限性所带来的问题。导出了变截面梁单元的单元刚度矩阵、单元等效节点荷载和单元一致质量矩阵。在此基础上,利用编制的程序进行了算例验证与分析。算例验证了本文理论的正确性,表明本文方法具有很高的计算精度。
1 引言
变截面梁在土木、航天和机械等领域有着广泛的应用,其优点是可以随内力的变化合理地进行截面设计,节省材料同时又减轻自重。因此,对变截面梁在荷载作用下反应数值求解方法的研究具有重要的理论意义和现实的工程背景。在诸多分析变截面梁反应的方法中,有限元方法是最具灵活性和适用性的方法之一。然而以往使用有限元方法进行的研究采用的是基于单元位移插值函数建立有限元公式的方法[1-3],由于位移插值函数不能有 效地描述单元内部的位移变化,所以此类方法对于变截面梁的分析存在着本质的缺陷。文献[4,5]首先求出截面变化具有简单规律(线性、抛物线等)的梁内部位移的解析解,然后使用将解析解作为单元位移插值函数的方法取得了好的计算结果,然而这类方法仅适用于梁截面的变化极其简单的情况。近年来,文献[6-8]中针对等截面梁提出了一种使用力插值函数的梁单元模型,并在解决材料非线性问题上显示出了巨大的优越性。该方法有限元公式的建立是基于单元内部力与单元节点力满足平衡。因该方法可以回避直接对单元内部位移场进行描述的问题,所以为更好地解决变截面梁的反应分析问题提供了新的思路。一些研究者[9,10]已采用这种思路对变截面梁进行了初步的研究。此外,另有一些研究者[11,12]在对变截面梁进行研究时虽未直接说明,推导公式时也是在基于首先使单元平衡条件得到满足这一前提下进行的。但这些研究多为静力分析,涉及动力分析的较少,并且在单元等效节点荷载和单元质量矩阵的构造方法上存在着一些问题。
本文根据文献[6-8]中的理论推导了变截面梁的有限元公式,并且给出了单元刚度矩阵、单元等效节点荷载和单元一致质量矩阵的建立方法。在此基础上,利用编制的计算机程序进行了算例验证和分析。
2 位移插值函数的局限性
传统梁单元在公式建立时使用事先给定且形状固定不变的位移插值函数来确定单元内部的位移变化。通常采用一次多项式构造轴向位移场,Hermitian三次多项式构造横向位移场。由于位移插值函数与几何相容方程结合可以得到单元内部任一点的应变与单元节点变形之间的关系,所以此种方法建立有限元公式的实质是首先使单元的几何相容条件得到满足。考虑平面弹性梁单元且仅受节点荷载作用的情况,对应等截面和变截面时分别有式(1)和式(2)成立:
式中qx和qy为梁轴向和横向分布荷载,E为弹性模量,I和I(x)分别为等截面梁和变截面梁的截面惯性矩,A和A(x)分别为等截面梁和变截面梁的截面面积,u(x)和v(x)分别为梁在x处的轴向和横向变形。
式(1a,b)的解分别为一次多项式和三次多项式,即单元内部位移场可以用一次多项式和三次多项式精确描述,所以传统梁单元使用一次多项式和Hermitian三次多项式构造位移场是精确的。然而当梁具有变截面时,式(2)表明一次多项式和三次多项式已经不能精确描述位移场,所以此时传统方法是近似的。由这种近似产生的误差为有限元的离散误差,表现为将实际物体离散成若干个单元时单元内部位移变化与实际物体位移变化不符,而将实际物体内部位移变化模式在进行有限元离散后代之以强制性的位移插值函数。在使用传统方法推导有限元公式过程中,位移插值函数发挥着至关重要的作用,最终影响单元刚度矩阵、单元等效节点荷载和单元一致质量矩阵的形成。所以,位移插值函数在描述单元内部位移场时的不准确性使传统有限元方法在分析变截面梁时表现出本质的局限性。在实际分析变截面梁时,经常需要采用细分单元的方法来尽量减少这种误差。此外,对于不使用一次多项式和Hermitian多项式构造位移场的梁单元,如使用自适应有限元方法[13],只是对单元内部位移场的进一步近似,这种离散误差仍不能被消除。
3 直接基于平衡的梁单元公式
3·1 单元平衡条件
单元平衡矩阵与单元截面面积的变化无关,所以不受截面变化的影响。将实际物体进行离散,单元内部的力场如用单元平衡矩阵构造则可以完全消除传统梁单元用于变截面梁时产生的离散误差。用这种方法建立有限元公式的实质是事先使单元的平衡条件得到满足,所以可以称之为直接基于单元平衡的方法。根据梁的微分方程、边界条件和虚力原理,可以得到x处截面力增量与单元广义节点力增量之间的关系:
积和惯性矩是x的函数,而不是常数。
式(8)中的K是3×3阶的,矩阵集装时需转换到整体坐标系中。考虑到广义坐标系中单元节点力在虚单元节点位移上所作的虚功等于整体坐标系中单元节点力在虚单元节点位移上所作的虚功,可以得到将K转换到整体坐标系中的转换矩阵,进一步得到整体坐标系中的单元刚度矩阵:KG=T(θ)TKT(θ) (9)式中KG为整体坐标系中的单元刚度矩阵。转换矩阵T(θ)见下式,其中s=sin(θ),c=cos(θ),θ为单元局部坐标系与整体坐标系之间的夹角,规定逆时针方向为正。
当梁为等截面时,用式(19)可得到与用传统方法相同的单元一致质量矩阵。对于整体坐标系中的单元质量矩阵,使用常用的由单元局部坐标系到整体坐标系的转换矩阵即可。值得一提的是单元等效节点荷载也可以用Φ与单元内部荷载q相乘然后沿单元积分的方法求出,与式(13)得到的结果是一致的。另外,文献[9,10]也曾提出过求单元一致质量矩阵的方法,思路均是先求单元内部位移场变化,然而使用文献[9]中的方法得到的单元一致质量矩阵与单元坐标系起点的选择(如图1中的I或J)有关,文献[10]由于选择了两端固结的梁作为虚力系统使其方法的推导十分复杂。本文给出的方法力学意义明确,避免了以上方法的不足。
4 算例分析
一竖直悬臂梁,长L=10 m,宽b=0·5 m,底端截面高h0=1 m,梁高以h(y) = h0-2y/25变化,材料弹性模量质量密度在静力分析中考虑受荷情况为垂直梁轴线的分布荷载q =200 kN/m。计算时采用三种方法进行有限元离散———本文方法的变截面梁、Ansys中的Beam54单元(此单元属传统梁单元,并可同时模拟变截面梁和等截面梁)。使用等截面Beam54单元进行模拟时,单元性质取整个单元两端点性质的平均值。图2给出了结构顶点横向位移和转动的计算结果,图3给出了梁的第一阶自由振动频率的计算结果。
需指出的是,由于算例中梁的截面变化较剧烈,用1~3个Ansys变截面梁单元进行模拟时,Ansys软件无法得到分析结果,所以图中对应于Ansys变截面梁单元的结果是从4个单元开始的。顶点位移和转动的解析结果分别为使用本文方法1个单元即可获得达到此精度的结果,而Ansys变截面单元和等截面单元分别需要超过110个单元和130个才能达到和解析解相同的精度。对于第一阶频率,使用本文方法4个单元得到的频率为10·24Hz,而使用Ansys变截面单元达到相同的精度需要30个单元。Ansys等截面单元对频率的分析精度不好,使用200个单元时仅能收敛到10·22 Hz。如果对求解精度要求不高,使用的Ansys变截面或等截面单元的数量可以大幅度减少,这一点可以从图2和图3中看出,但同时也可看出本文方法仍是具有非常大的优势。
5 结 论
本文采用直接基于单元平衡的思想给出了计算变截面梁静力及动力反应的方法。使用直接基于单元平衡的思想推导有限元公式,不依赖于单元位移插值函数,避免了单元位移插值函数用于近似变截面梁内部位移变化时不准确的局限性。
给出了单元刚度矩阵、单元等效节点荷载和单元一致质量矩阵的推导方法。算例验证了本文理论的正确性,表明本文方法可显著提高对位移反应和频率的求解精度。由于使用较少数量的单元即可获得较高的精度,因此也可使分析问题的计算效率得以提高。
