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

1 概述

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

（1）   基础理论

（2）   商软操作

（3）   自编程序

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

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

                                            

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

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

        

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

==第40篇：Abaqus、iSolver与Nastran梁单元差异（4）形心、剪心和偏置==

从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梁差异（4）-形心、剪心和偏置。

                                             

1.1 形心、剪心和节点坐标

一根梁虽然在有限元内部计算时采用线单元表示，但实际上只是把梁的三维形状简化到了参数中，如果显示一个梁的具体形状，譬如典型的T型梁：

将会有下面三个特殊点：

（1）形心C

就是的二维截面几何图形的重心。

（2）剪心S

剪心为两个截面方向的集中力都不会导致扭转的交点。譬如槽钢的剪心S大体如下，显然和形心C不是同一点：

如果在y方向加集中力，可发现只有通过剪心才能不发生扭转，即下方的第二图：

 

（3）节点位置N

节点位置就是建模时输入的位置，譬如下面的三个节点就是0,0,0、100,0,0,和200,0,0

那么在CAE软件中就是线显示的位置：

 

• Nastran中节点称为Grid，对应bdf中的GRID关键字

既然梁简化为线单元的过程中有三个特殊点：形心C、剪心S和节点N，那么问题来了：

（1）   输入的力加载在什么位置？

（2）   计算采用的运动的相对位置是什么？

1.2 输入的集中力和力矩加在哪个点

梁单元截面上的加载就是集中力、yz方向的弯矩，x轴的扭矩，我们分开讨论：

1.2.1 集中力

对于集中力，实际生活告诉我们一根手指去压迫梁截面，那么手指和截面接触的位置显然就是加载点。

实际力和力矩可以加载到任意点上，只要和实际的力和位移关系一致就行，但有限元让用户输入集中力和力矩时没有输入截面坐标，只可能默认是上面三个特殊点之一。

具体加载到什么位置可以做个简单的实验证明：

Example1：Ex1-LoadPosition\

譬如上面的T型材，如果将节点位置N改到剪心S上方（实际不是改节点位置，而是改偏置）

将左端固定，右端在z方向加集中力载荷。

那么，可能如下：

（1）   如加载到形心C，右侧节点2绕x轴逆时针旋转，得到扭转角度UR1<0

（2）   如加载到节点位置N，右侧节点2绕x轴顺时针旋转，得到扭转角度UR1>0

（3）   如加载到剪心位置N，不旋转，得到扭转角度UR1=0

在iSolver中计算如下：

可发现UR1>0，说明是加载到节点上的。

1.2.2 力矩

力矩在有限元中和力一样都是加载到过节点N的转动轴上，注意力矩是输入，与梁实际的转动轴无关。譬如下图节点如果在底部，那么力矩M通过力F*到节点的距离应该表示为下图所示：

• 但我们找不出像力一样可以证明力矩加载位置的例子，哪位读者要能证明一样也可以联系我们。

我们查看Nastran的帮助手册，也可发现Nastran的集中力也是加载在节点Grid上。

所以，集中力和力矩就是加载在节点Node所对应的N（x、y、z）上，只不过N（x、y、z）可能恰好是剪心，也可以是形心，当然，也可能是不在这两个点上的任意值。

1.3 计算采用的运动的相对位置

力的加载位置和梁的真正后台计算的相对位置不同，梁的运动分为平动、绕x轴的扭转、绕y/z轴的弯曲转动，平动显然是节点本身的运行，不做讨论，但转动是绕转动轴运动的，这个转动轴和输入的力矩认为转动轴不是一回事，由上面讨论输入的力矩认为转动轴总是过节点N的，但这边运动的转动轴是软件内部计算采用的转动轴，每个软件还是不同的。

1.3.1 扭转的扭矩轴

有限元中，如果受到剪心位置的扭矩，那么扭转角沿过剪心的x方向轴计算。

即J中的距离是相对剪切中心来计算。

当剪心不在节点上时，无论材料力学理论还是Abaqus/iSolver/Nastran软件会认为扭转轴还是绕过剪心的轴线，而不是过节点的：

此时，显然J不变，但需要把节点的扭矩和扭转角变换到剪心。

当一个集中力矩移动到另一点时，力矩和扭转角度都不会变，但绕剪心的旋转对节点导致截面两个方向的位移v、w的变换，且显然与剪心和节点的距离NS有关，具体这个v、w的变换如何影响刚度也可看我们以前梁单元刚度阵的文章。

 

1.3.2 弯曲的转动轴

弯矩是绕某个轴的旋转，我们来分析究竟绕过节点还是形心旋转？

1.3.2.1 一般的材料力学理论的推导过程

我们先从一般的弯矩公式的最原始的推导说起，在受弯矩作用下，梁做纯弯曲变形，整个梁的每条平行与轴线方向的纤维都弯曲为圆形。

那么纤维SS’的相对nn’的升长量是

如果nn’是中性轴，也就是nn’的长度就是纤维SS’原来的长度，那么就是SS’的应变。由应力应变关系：

内外力矩平衡：

其中：Iz为

得到应力：

材料力学书中一般认为y表示所在截面点到中性轴的y位置，也就是说弯曲是绕中性轴的运动。

但从上面的推导可以看出，要上面的公式成立：需满足：

1.纯弯曲情况。

2.满足中性轴nn’是纤维SS’的没有载荷时的原始长度。

3.弯曲后截面依然垂直中性面，否则下面等式不成立，也就是说只适合Euler梁。

1.3.2.2 有限元中的力矩弯曲轴

实际的情况很多时候都不是纯弯曲，同时中性轴也会拉伸，所以到底绕哪个轴弯曲每个软件也不同，Abaqus中取的是过节点的弯曲轴，取过节点的好处是由于上面说到的集中力和位移等都是节点的物理量，这样所有的物理量都是节点上的，不用变换。相当于截面上力的内矩采用下图表示：

而Nastran中弯曲轴就不同了，取的是形心的弯曲轴。类似于软件内部要做两步操作：

（1）   先将计算转动轴放到中性面：

 

（2）和上面的扭转一样，Nastran的梁单元在后台需要将形心的位移转换到节点上，譬如下方在形心处的弯曲角度θ会导致x方向的位移:

1.4 实际工程上的弯曲轴

实际的弯矩外力是加在过节点Node还是形心的轴向的呢？也就是说梁的截面是沿哪个轴弯曲的呢？个人没做过实际的物理试验，只是查了一下力矩的加载方法，一种加载方式是在梁的两个边缘加力，然后用力乘以离弯曲轴的距离当做力矩。但试验室截面到底是观测到过节点还是过形心个人理解很难说的清楚，如果问一个小朋友，看上面的弯曲图，转动角度θ是绕哪个轴旋转，也许很多人就会说是底部，而不是中性轴，也就是说从外形上很难确定哪个轴旋转，或者你可以认为是绕任意一个轴旋转。

那么有限元中这两种方式模拟结果是不是一致呢？ 我们只以一个简单的例子来说明一下。

1.4.1 算例2

Example2：Ex2-TBeamNodeInShearCenter\

我们还是以T型材的一个梁单元为例。

先在iSolver中建模，为减少横向剪切刚度的影响，我们梁长度L取为240，仅取一个单元。

材料杨氏模量=1.5e6，泊松比=0.3。

 

节点在剪心处，此时T型材的截面尺寸和偏置如下：

1.4.2 工况和实测结果

右端全约束，取三种工况：

工况（1）：左端自由加力，加y-方向集中力F=1：

iSolver导出到Abaqus B31单元和Nastran的CBEAM单元绕z轴转动角度UR3结果都是5.02e6，如下：

 

• 上面说了弯曲的材料力学理论只适合Euler梁，但我们Abaqus还是采用了B31，仅仅因为B31更常用，但就算是采用B33，Abaqus内部计算依然采用过节点的转动轴。

工况（2）：左端自由加力矩，加z+方向力矩M=1：

iSolver导出到Abaqus B31单元和Nastran的CBEAM单元UR3结果都是2.09，如下：

 

工况（3）：在工况（2）的基础上左端简支

Abaqus B31单元和Nastran的CBEAM单元UR3结果分别是5.229和4.24，如下：

1.4.3 结果分析

明显工况（1）和（2）采用经过节点的弯曲轴的Abaqus和采用经过形心的弯曲轴的Nastran结果基本一致，但工况（3）明显差异较大。我们没有严格的证明，我们猜测无论采用哪个弯曲轴都是正确的，只是多一个平动位移的变换，都是严格满足梁的外力矩=内力矩，外力=内力的关系的。所以（1）和（2）的结果一致，但如果像工况（3）把平动位移约束，弯曲轴的不同导致的差异可能就非常大了。

进一步，对于工况（3），Abaqus和Nastran那么哪个更正确呢？或者谁更可信呢？当我们划分为100个梁单元时，发现Abaqus和Nastran的UR3结果均是4.2左右，我见过这个问题的理论值都是建立在节点在形心位置的，至于现在节点不在形心的情况，没找到对应理论值，只能猜测4.2就是比较精确的实际工程值，那么回头看如果仅仅是一个单元，那么Nastran更精确，当单元足够多时，Nastran和Abaqus结果都一致。

 

• 对于Abaqus的B31一个单元的结果不如Nastran的情况我们发现过很多次了，譬如简单的悬臂梁加弯曲力的例子，猜测Nastran的梁理论是精确解，无论单元个数多少。但Abaqus是近似解，只有单元多了才和理论解一致，这两种情况工程上都可用。我们没有找到对应的理论解释，如果哪位高手知道为何这样，欢迎联系我们。

1.5 偏置

最后我们聊一下形心和剪心位置相关的偏置。

形心和剪心具体是在截面几何图像上哪个点完全由几何形状决定，但形心点和剪心点在三维空间的坐标与这个梁的空间摆放位置有关，在有限元中，这个空间摆放位置由两个参数决定：

（1）   节点N的位置决定整体几何截面的2D面。

（2）   偏置Offset决定形心/剪心相对节点N的位置。

节点N相对简单，用户设置多少就是多少，但偏置Offset每个软件中的设置不同。

1.5.1 Abaqus

Abaqus的偏置只有一处地方，就是在设置型材的时候设置Offset的参数，譬如T型梁中设置I：

 

1.5.2 Nastran

Nastran的Offset有两处：

（1）   在单元Property时设置Offset，它默认是全局坐标系的，后台是对CBEAM和CBAR关键词修改，表示单元的偏移。

（2）单独在PBEAM中设置Offset，此时Offset是相对型材截面局部坐标系的，表示形心的移动，在BDF中PBEAM关键词中可以修改，没找到界面的修改地方，个人基本不用，不是很懂。

1.5.3 iSolver

Nastran相比Abaqus的优势：

（1）Nastran所有的型材都可以设置偏置，而Abaqus很多型材如果要偏置就非常麻烦，譬如Rectangle没有偏置，只能采用Trapzoid中设置上下底相等或者采用Arbitrary梁。

（2）看上面解释可对形心和节点分别设置偏置，个人想不出应用场景，总觉得型材安装位置应该只能整体移动，而不是节点和形心可以分别移动。

Nastran相比Abaqus的劣势：

很多实际情况是型材和下面的安装的蒙皮的相对位置是定的，但在三维空间的位移是变化的，所以个人觉得Abaqus的相对坐标系的偏置更实用一点。

也正是基于Abaqus和Nastran的优缺点考虑，iSolver采用了Abaqus的局部坐标系的偏置，同时和Nastran一样可以对矩形或者L型等设置双向偏置。当然，更好的设置偏置的方式是自动偏置，这是Nastran和Abaqus均没有的功能。

1.6 视频讲解和操作验证演示

如果觉得上面的文字太复杂，也可以看一下视频的简要讲解和软件演示，地址如下：

https://www.jishulink.com/college/video/c12884 20理论系列文章40-梁单元差异(4)-形心、剪心和偏置

==以往的系列文章==

 ========第一阶段========

第一篇：S4壳单元刚度矩阵研究。

http://www.jishulink.com/content/post/338859

第二篇：S4壳单元质量矩阵研究。

http://www.jishulink.com/content/post/343905

第三篇：S4壳单元的剪切自锁和沙漏控制。

http://www.jishulink.com/content/post/350865

第四篇：非线性问题的求解。

http://www.jishulink.com/content/post/360565

第五篇：单元正确性验证。

https://www.jishulink.com/content/post/373743

第六篇：General梁单元的刚度矩阵。

https://www.jishulink.com/content/post/403932

第七篇：C3D8六面体单元的刚度矩阵。

https://www.jishulink.com/content/post/430177

第八篇：UMAT用户子程序开发步骤。

https://www.jishulink.com/content/post/432848

第九篇：编写线性UMAT Step By Step。

http://www.jishulink.com/content/post/440874

第十篇：耦合约束（Coupling constraints）的研究。

https://www.jishulink.com/content/post/531029

========第二阶段========

第十一篇：自主CAE开发实战经验第一阶段总结。

http://www.jishulink.com/content/post/532475

第十二篇：几何梁单元的刚度矩阵。

http://www.jishulink.com/content/post/534362

第十三篇：显式和隐式的区别。

http://www.jishulink.com/content/post/537154

第十四篇：壳的应力方向。

https://www.jishulink.com/content/post/1189260

第十五篇：壳的剪切应力。

https://www.jishulink.com/content/post/1191641

第十六篇：Part、Instance与Assembly。

https://www.jishulink.com/content/post/1195061

第十七篇：几何非线性的物理含义。

https://www.jishulink.com/content/post/1198459

第十八篇：几何非线性的应变。

https://www.jishulink.com/content/post/1201375

第十九篇：Abaqus几何非线性的设置和后台。

http://www.jishulink.com/content/post/1203064

第二十篇：UEL用户子程序开发步骤。

https://www.jishulink.com/content/post/1204261

========第三阶段========

第二十一篇：自主CAE开发实战经验第二阶段总结。

https://www.jishulink.com/content/post/1204970

第二十二篇：几何非线性的刚度矩阵求解。

http://www.jishulink.com/content/post/1254435

第二十三篇：编写简单面内拉伸问题UEL Step By  Step

http://www.jishulink.com/content/post/1256835

第二十四篇：显式求解Step By Step。

https://www.jishulink.com/content/post/1261165

第二十五篇：显式分析的稳定时间增量。

http://www.jishulink.com/content/post/1263601

第二十六篇：编写线性VUMAT Step By Step。

https://www.jishulink.com/content/post/1266640

第二十七篇：Abaqus内部计算和显示的应变。

https://www.jishulink.com/content/post/1273788

第二十八篇：几何非线性的T.L.和U.L.描述方法

https://www.jishulink.com/content/post/1282956

第二十九篇：几何非线性的T.L.和U.L.转换关系

https://www.jishulink.com/content/post/1286065

第三十篇：谐响应分析原理

https://www.jishulink.com/content/post/1290151

========第四阶段========

第三十一篇：自主CAE开发实战经验第三阶段总结

https://www.jishulink.com/content/post/1294824

第三十二篇：谐响应分析算法

https://www.jishulink.com/content/post/1299983

第三十三篇：线性瞬态动力学

https://www.jishulink.com/content/post/1302074

第三十四篇：非线性瞬态分析

https://www.jishulink.com/content/post/1787283

第三十五篇：接触求解算法

https://www.jishulink.com/content/post/1792869

 第三十六篇：DLOAD用户子程序开发步骤

https://www.jishulink.com/content/post/1826803

 第三十七篇：梁单元差异（1）-理论基础

https://jishulink.com/content/post/1872208

 第三十八篇：梁单元差异（2）-梁截面方向

https://www.jishulink.com/content/post/1874628

 第三十九篇：梁单元差异（3）-剪力和弯矩

https://www.jishulink.com/content/post/1876013