有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础

SnowWave02 浏览：2862 评论：20 收藏：21

（原创，欢迎转载，转载请说明出处）

1 概述

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过

（1）基础理论

（2）商软操作

（3）自编程序

三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。

有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图3

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元求解器，通过自研求解器和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。

通用结构有限元软件iSolver介绍视频：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图5

http://www.jishulink.com/college/video/c12884

==第37篇：Abaqus、iSolver与Nastran梁单元差异（1）-理论基础==

从NASA在1964年推出Nastran算起，结构CAE软件已经发展了将近60年，虽然在60年内商用结构软件层出不穷，但最流行的通用结构CAE软件现在依然只有Nastran、Ansys和Abaqus三款软件，iSolver以前的结果精度完全对标Abaqus，线性和材料非线性与Abaqus计算结果对比精度已经没有误差，但在推广过程中，发现在线性问题上，客户只相信Nastran的结果，因为很多工程的校核规范和经验方法都以Nastran的结果为基准，如果换成其它软件，那么那些经验系数和校核公式都可能会修改，而这些修改还需要做大量的工程验证，这是客户承受不了的代价。客户改不了的习惯只能是自主软件改，我们的自主软件只是一个后来者，想要推广只能按客户要求修改精度，因此，iSolver开始兼容Abaqus和Nastran的精度，计算出的结果既可以对标Abaqus，修改设置后又可以对标Nastran，因此需要对Abaqus和Nastran的各方面的算法差异进行深入研究。如果是线性问题，那么Nastran和Abaqus的精度误差主要体现在单元算法、边界处理、MPC约束关系等，在2017年第二篇：S4壳单元质量矩阵研究文章中我们就曾经分析过Abaqus的S4壳单元和Nastran的Quad4壳单元质量矩阵的内部实现方式和差异，在这里主要研究Abaqus、iSolver与Nastran梁单元差异，由于这三款软件的梁单元的差异较多，我们分几篇文章来说明，本篇是Abaqus、iSolver和Nastran梁差异（1）-理论基础。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图8

1.1 梁单元基础理论

1.1.1 有限元中梁单元基本假设

在有限元中，梁单元有个基本假设，就是原来的截面（譬如下面的红色表示的截面）受弯曲力后依然是平面。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图10

那这个假设和实际符合吗？

如果是在纯弯矩情况下，由试验得到弯曲时依然保持平面，符合实际：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图12

如果一个橡皮条做成的悬臂梁，那么原来在铅直线上画的直线受力后将会变成如下图所示的m1n1的曲线。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图14

很明显，上面的弯曲下原来的平面mn已经变为了曲面m1n1了。如果是这样，不就和梁弯曲后依然为平面的假设矛盾了吗？

这个就得分清平面假设的目的了，这个假设的目的是得到弯曲时正应力的分析情况，即

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图16

正应力随y方向的变化为：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图18

很多工程试验表明，只要按上面公式求解，就能得到和实际工程问题符合的位移/应力等结果，因此问题不在于截面弯曲后是否还是平面，而在于弯曲后是否还能用上面的公式。

上面悬臂梁受力的例子，截面是已经不再是平面了，如下图：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图20

但此时可以发现所有横截面的翘曲都一样，导致m1p1=mp，即纵向纤维上产生的伸长或缩短不受剪力影响，那么由纯弯曲导出的正应力公式依然成立。

1.1.2 梁单元理论分类

在弯曲时，梁的下端面变短被压缩了，上端面变长被拉长了，那么中间必然有个长度不变的平面，就称为中性面。受力前截面和中性面是垂直的，那么弯曲后的截面和中性面就有两种可能：

（1）如果是细长梁，那么可以认为依然垂直，此时就称为Euler梁。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图22

（2）如果是厚梁，那么不再垂直，此时就称为Timoshenko梁。显然Timoshenko梁覆盖了Euler梁的功能，在细长梁的情况应该就和Euler梁完全一致。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图24

分类是极其简单的，但也许有人会认为Euler梁就等价于三次梁，即用三次形函数来表示位移，而Timoshenko梁是一次梁，即用一次线性函数来表示位移，个人理解不是这样的。

（1）首先是梁包括轴向拉伸、扭转、横向弯曲、剪切行为等，有限元中对不同行为都是用不同次数的形函数的，上面所说的一次和三次仅仅是横向弯曲时采用几次多项式的问题。

（2）无论是Euler还是Timoshenko梁，形函数不一定非要用一次和三次，你完全可以用三次形函数来表示Euler梁，或者用五次形函数来表示Timoshenko梁，甚至压根不是用多项式直接来做刚度阵中的形函数，譬如很多人算变截面梁的刚度阵时采用的是基于力的形函数，先计算柔度矩阵，再计算刚度矩阵。

计算横向剪切应变时：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图26

Euler梁和Timoshenko梁区分的唯一标识是有限元中是否考虑横向剪切效应。如采用Euler梁，那么显然

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图27

使得横向剪切应变

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图30 也就没有剪切效应。而Timoshenko计算出来的横向剪切效应就不能忽略了。

1.2 Abaqus的梁单元

Abaqus最常用的梁单元类型为B33和B31，分别对应Euler和Timoshenko梁。在Mesh模块的ElementType（就是有个SAR图标的按钮）中设置，如下。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图32

不过Abaqus界面上的看的有点让人糊涂，shear-flexible是Timoshenko梁，Cubic formulation是Euler梁，要么统一用是否可以剪切，譬如shear-flexible和no shear-flexible，要么统一用弯曲形函数的多项式次数，譬如linear formulation和cubic formulation，也可能是我没明白Abaqus的用意，要有大神了解的话也说一下。

1.3 Nastran的梁单元

1.3.1 Nastran梁单元的设置

Nastran的梁单元主要有CBAR和CBEAM两种。

（1）在Patran界面上，在Property模块下设置梁单元属性，默认情况下Option为General Section/Standard Formulation：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图34

输出到Nastran的Bdf中，可发现单元类型为CBAR：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图36

（2）在Patran界面上，在Property模块下梁单元属性Option改为General Section（CBEAM）/Standard Formulation：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图38

输出到Nastran的Bdf中，可发现单元类型为CBEAM：

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图40

我以前想当然的认为CBAR和CBEAM分别就是Euler梁和Timoshenko梁，后来深入了解后，发现根本不是这回事。CBAR和CBEAM即是 Timoshenko梁，也是Euler梁。

部分Nastran理论说明介绍，默认情况下CBAR没有设置K1，K2，得到的是Euler-Bonourli梁，而CBEAM如果设置了K1，K2，是Timoshenko梁，但我们实测发现其实他们都默认设置了剪切效应。

在Patran梁截面设置上，有两个的参数Shear Stiff.Y和Shear Stiff.Z，这两个就是通常认为的横向剪切效应因子K1，K2。对于矩形它的默认值就是5/6。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图41

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图43

Abaqus也有K1、K2，但值和Nastran不一样，譬如矩形Abaqus取的是0.85， 5/6=0.8333是有严格的理论支撑的，正常来说这才是正确解，Abaqus为何取0.85不是很理解。如果你编的自主软件结果想要对标哪个商业软件，就只能按这个商软的参数修正，当然像iSolver精度要对标Nastran和Abaqus两款软件，后台就得判断取不同的K1、K2。
Patran界面上可以查看几何参数组成的截面最后计算时的面积、形心位置、弯曲惯性矩、扭转惯性矩具体结果，实测发现Nastran和Patran计算出来的面积、形心位置、弯曲惯性矩、扭转惯性矩并不完全一致。

CBAR和CBEAM如果设置了K1，K2，那么考虑了剪切，为Timoshenko梁，如果没有设置K1，K2，得到的是Euler-Bonourli梁，Patran中默认情况下设置了K1、K2，所以默认情况下都是Timoshenko梁。同时，无论是哪种梁单元，Nastran的结果在细长梁下都会和Euler梁非常接近。

1.3.2 Nastran的梁单元区别

既然CBAR和CBEAM都是Timoshenko梁，那他们的区别是什么呢？我们的理解，最大的区别是CBAR只能支持剪心和形心相同的情况，CBEAM支持剪心和形心不相同的情况。

还有其它的区别是CBAR只能常截面梁，CBEAM能处理变截面梁等，但这不是最主要的区别。

对剪心和形心相同的情况，那么用哪种单元都没问题。但对不在同一点的情况，如果采用CBAR就可能导致结果精度有误差，譬如下面的四种图形，显然前三种剪心和形心明显不同，对扭转、弯矩、不同方向的集中力等复杂加载，此时如果采用CBAR就会导致没有剪心和形心不同心导致的耦合效应，结果不一定正确。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图45

Patran在帮助文档中已经明确说明了该问题，同时对L型材，产生的梁单元类型都自动设为了CBEAM，但对T型材还是可以采用CBAR，这也导致很多采用CBAR来模拟T型材的结果和实际有较大偏差。

1.4 iSolver的梁单元

iSolver也包括Timoshenko和Euler两种梁单元。界面上设置就完全按照是否包括剪切效应来决定后台采用Timoshenko还是Euler梁。

有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图47

其中：

（1）如勾选Include transverse shear，那么是Timoshenko梁，梁单元完全按照Abaqus的理论实现，有兴趣可看有限元理论基础及Abaqus内部实现方式研究系列6：有限元理论基础及Abaqus内部实现方式研究系列37：梁单元差异（1）-理论基础的图48 General梁单元的刚度矩阵的文章。

（2）如不勾选Include transverse shear，对Euler梁，一开始我们也准备按Abaqus实现，但后来深入研究后发现Abaqus的B33虽然是Euler梁，它的后台算法和普通的3次梁还是做了较大的修正，它的刚度阵是14X14，比普通的12X12多了两个内部自由度，不清楚为何这么做，没找到它对应的理论，同时，我们对比发现当梁数目不多时，Abaqus的B33精度并不很好，如下面用户采用Abaqus/Nastran/Ansys/Midas/iSolver五款CAE软件计算的球面网壳结构的自振频率的例子，当两个结点间采用一个欧拉梁单元时，可发现Ansys、Nastran、iSolver、Midas四者都很接近，均为10.0x，差异仅为0.5%内，但Abaqus的B33却得到9.4386，和前四者都差6%左右，因此iSolver的Euler梁没有按照Abaqus的精度来实现，而是按照Nastran的算法实现，导致iSolver的Euler梁和Nastran精度完全一致。