有限元理论基础及Abaqus内部实现方式研究系列45：约束关系（1）-统一形式
（原创，转载请注明出处）
1 概述
本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过
（1） 基础理论
（2） 商软操作
（3） 自编程序
三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。
有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。
一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元软件iSolver，通过自研CAE软件和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。iSolver包括完整的前后处理和有限元求解器，功能如下，有兴趣可直接在下面网址下载：
百度网盘链接: https://pan.baidu.com/s/10d6jHdZ01SBY2JxiS6bffw 提取码: 6fdf
2 约束关系
我们在系列文章第35章介绍了接触分析及常用的基本算法。现在Abaqus、LS-DYNA、Ansys等结构商软都说可以处理复杂的上万零部件接触的整车、整机等模型仿真，没做过实际的这种仿真分析，很好奇，接触分析算法往往涉及大变形、边界不连续，只要输入条件或者算法稍微变化一些，两个零部件算出来的接触结果就可能差异很大，更不用说上万个零部件的接触结果了，对这种大规模组装模型的仿真结果不知如何来判断它的可靠性，像普通的只校核一下材料的应力还是看一下动画是否和试验一致？毕竟仿真只有简单的标准来判断结果的正确性才能在企业中起到真正辅助设计的作用，如果你恰好做过，不妨也简单介绍一下你的经验。对自研CAE软件开发者来说，自研结构CAE软件是否真的要和商软去比拼接触等复杂算法还是花更多时间在精雕细琢那些常用功能上，这也是开发者需要慎重考虑的问题，而且很多自主CAE软件连常规线性问题都算的不对，或者都没法用鼠标稳定的走完那些材料、属性、边界、加载等流程，用户又怎么会相信你能算对接触这种复杂问题的？
不管怎么样，从有限元实现的角度来讲，如果想做真正实际工程中的接触分析，那么首先需要去研究约束关系，接触分析在有限元中也仅是约束关系的一种。有限元中的约束很多场景大家用的是边界中的简支、固支等约束，但从更广泛的角度上讲，只要表示一个节点的某个自由度依赖于其它的节点自由度或者取某个特定值，就可以称为约束关系。只不过对固支、简支等直接自由度=0，在有限元中直接减缩刚度阵就行，很容易求，但对节点自由度相互依赖的约束关系就比较复杂了。约束关系主要有两类。
（1） 一类是MPC点之间的约束。Nastran的MPC的灵活度要远远超过Abaqus，Nastran的主节点可以选择123自由度，也可以对每个从节点设置不同的自由度，还能主节点和从节点互相包含，Abaqus更多的是只负责80%的常用应用场景，复杂功能让你编子程序，但事实上一线仿真工程师又有多少人愿意编子程序呢？这种做法导致虽然Abaqus无论从用户体验、非线性还是商业化都比Nastran好很多，但很多线性的工程复杂问题还是没法替代Nastran。
（2） 另一类是Contact、Tie等的面之间的约束关系。在这方面Abaqus要明显强于Nastran了。
我们将用统一的公式来求解这两类关系，同时也从软件实现层面说明一下针对这两类情况的各自差异。分几篇文章来介绍约束关系，本篇是约束关系（1）-统一形式，既然接触仅是约束关系的一种，那么MPC、Tie、接触等的求解过程也是很类似的，这里将介绍一下这些约束关系如何表达为统一形式。
3 统一形式的约束关系
在没有约束关系时，如下图情况，物体在体外力和面外力作用下变化。
有限元方程按照虚功原理求解，在物理上可解释能量守恒原理，即在某一个时刻点，假定在外力作用下有个虚拟的位移，那么外力在虚拟位移下做的虚功=内部应变能的变化相同。
虚功原理中的每项都表示各自区域在虚位移下的能量变换
（1） fv是每单位体积内的力，外力fv和位移相乘表示单位体积内的虚功，所以对体积积分
（2） fs是每单位面积上的力，外力fv和位移相乘表示单位面积上的虚功，所以对面积积分，推论就是最后一项应力和应变是单位体积内的内能，所以对体积积分。
当存在约束关系时，在能量中加入约束关系相关的一项：
显然的单位也是该约束关系所在区域在虚位移下的能量单位。
约束关系在有限元中可以分为两大类：
3.1 点之间的约束关系
（1）点之间约束关系，最常见的是节点之间的刚性连接，Nastran中称为RBE2，在Abaqus或者iSolver中称为Kinematic Coupling，此时可以认为Master节点和Slave节点之间焊死在了一个刚性无穷大的直杆上。在实际情况中，Slave节点之间没有相对位移，但由于计算时很多时候默认为小位移，反而导致Slave节点之间是有相对位移的。
此时：
其中，
为Master节点的旋转向量，为从Slave节点到Master节点的向量。
这种节点之间的约束关系还有一种常见情况就是节点间的分布耦合连接RBE3，将施加在Master节点上的力和力矩按照加权系数分配到Slave节点上，从而实现载荷在单元间的传递，此时就避免了RBE2太刚的问题，典型应用场景譬如模拟圆管对壁面的压力作用。
分布耦合连接是将Master节点的力和力矩按某种规则分配到Slave节点，保证力和力矩分别相等，即：
上述对六自由度的Master节点来说，一共只有六个方程，但每个Slave节点都有6个位置量，所以对多个Slave节点情况解不唯一，Abaqus、iSolver、Nastran等都有各自不同的分配原则，一般都是假定Slave节点不再存在Master节点分配过来的弯矩Ms，同时，Fs等比例分配。
对这种点的约束关系，虚功原理中增加的能量项表示为：
注意，因为是点之间的约束，所以不需要对面或者体积分，类似质量点对质量阵的贡献或者集中力对载荷向量的贡献时也不需要积分一样。
3.2 面之间的约束关系
面之间的约束关系，最常见的就是Master面和Slave面/点的接触连接关系。
在有限元软件中设置完接触关系对后，对Slave面上的任意一点xs，在Master面上寻找和xs距离最短的点，譬如xm点，此时xs和xm之间的距离为h，显然，Master面在xm点上的法向n指向xs点（要不然就不是最近邻点了）。也就是说slave面上任意一点在Master面上都有一个点来对应，这个点我们称为锚点，锚点的法向指向Slave点是一个必要条件（不是充分条件，因为有可能Master有两个锚点法向指向同一个xs点，譬如Master是个拐角的情况）。如下图所示。
此时距离(真实的距离肯定是正的，但我们这边为了方便取有向距离)
设接触面法向压力为p，对硬接触，当两个面之间的距离h是0时，压力p>0。有距离h<0时，没有接触关系，此时压力p为0。
即约束方程为
这是对面上的任意一点。计算机是无法处理面上所有点的问题的，所以也只能划分成有限单元，得到有限节点之间的关系。
划分完网格后，上述接触面上必然也是节点组成，如下：
对每个Slave节点xs，在Master接触面上寻找对应的锚点，因为Master面由不连续的节点组成，所以这个锚点很多情况都不在Master节点上。譬如上面图示，此时就是该点所在单元对锚点的插值，同样满足这个锚点的法向n(r,s)（也由所在单元的节点在r，s点插值得到）指向xs。
同样，此时的约束方程为：
这种面面之间的一种特殊情况就是Master面和Slave面/点的粘贴连接关系，表示Slave面/点用胶水粘在Master面上，所以在Nastran中称为Glue，在Abaqus中称为Tie，一般用于两个不同网格之间边界耦合在一起，譬如下方的桥身和桥墩，实际上Slave面/点是焊死在Master面的，Slave面/点会随着Master面一起拉伸压缩移动，但Slave面/点在运动过程中不会脱落，最终反应到刚度阵依然是固定的Slave节点自由度和Master节点自由度之间的约束关系。
对这种面与面/点的约束关系，增加的能量项表示为：
Sc为接触面。
上式和节点的约束关系相比，可知pms就是Lagrange因子。
在实际三维接触分析中，接触力除了法向n的压力还有两个切向v1和v2的摩擦力。得到实际接触力由三个正交的力分量组成：
简单起见，我们假定不滑动，此时由于接触关系增加的能量项为：
4 统一形式的约束关系
上述不同的仅是约束关系带来的积分范围和约束方程的个数和物理量。去掉所有的积分形式，我们可以统一写成如下形式:
其中上面的正负不影响结果，仅仅和h的正负有关。
5 更大统一形式的有限元的约束关系
MPC、接触等都是某部分节点和另一部分节点之间的关系，那么有限元中梁、壳、体普通的单元是什么呢？从更大的统一形式来说，这些普通单元也仅仅是节点之间的约束关系而已，只不过这个关系可能更加复杂点。
以两节点梁单元为例。
此时，每个节点6个自由度，那么单元刚度阵K为12X12的矩阵。K可以按6x6矩阵分块，每块小矩阵都比较类似，那么可以分成4块：
（1）左上角为第一个节点自由度相关的6X6的小矩阵做研究对象。iSolver软件中内部打印如下
K的任意一项Kij都是i和j两个自由度的乘积，其中对角元素为节点1和节点1自身的自由度的乘积，譬如K11就是u和u的乘积，而非对角元素都是耦合项，譬如K15就是u和θy的乘积。
（2）右下角类似，只不过是节点2和节点2自身的关系。
（3）左下角和右上角类似，是节点1和节点2的关系。
所以，可以总结有限元的所有单元，包括MPC、接触等最终计算的都是这个单元所涉及点两两之间的关系，这个关系的最终形式都体现在刚度阵上。有限元简单来讲就是求各个点之间关系的公式，只不过有些公式简单，譬如两个点之间用线性弹簧，有些复杂，譬如几何大变形时两个节点之间或者接触时两个面之间的关系。
6 以往的系列文章
6.1 ========第一阶段========
第一篇：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
6.2 ========第二阶段========
第十一篇：自主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
6.3 ========第三阶段========
第二十一篇：自主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
6.4 ========第四阶段========
第三十一篇：自主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
第四十篇：梁单元差异（4）-形心、剪心和偏置
https://www.jishulink.com/post/1888000
6.5 ========第五阶段========
第四十一篇：自主CAE开发实战经验第四阶段总结
https://www.jishulink.com/post/1905963
第四十二篇：声学分析(1)-有限元
https://www.jishulink.com/post/1912491
第四十三篇：声学分析(2)-边界元
https://www.jishulink.com/post/1923936
第四十四篇：声学分析(3)-湿模态
回复 
