有限元理论基础及Abaqus内部实现方式研究系列：约束关系-Lagrange因子法求解

1 概述

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

1） 基础理论

2） 商软操作

3） 自编程序

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

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

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

百度网盘链接: https://pan.baidu.com/s/10d6jHdZ01SBY2JxiS6bffw 提取码: 6fdf

2 约束关系

现在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了。

我们将用统一的公式来求解这两类关系，同时也从软件实现层面说明一下针对这两类情况的各自差异。分几篇文章来介绍约束关系，上一章我们介绍了基于点或者面之间的约束关系的形式是统一的关系，这章我们就将以此具体推导这种统一的约束关系在有限元中如何采用Lagrange因子法求解。

3 Lagrange因子法求解约束方程

没加约束前的虚功原理表示为：

4 各个软件的约束关系的具体实现方式

约束关系虽然原理都可以用Lagrange因子法来求，但软件的实现往往更复杂，实际编程有具体的两种实现方式：

1） 既然上一章说了约束关系可以类比普通的有限单元，那么可以把约束关系按单元来实现。

2） 约束关系也就是Master节点自由度和Slave节点自由度之间的系数，直接把这个约束关系写下来放到全局刚度阵就行。

4.1 Nastran的实现方式

Nastran理论手册中说明了一开始版本两者都有，但后面更多推荐用单元方式来实现，而且默认情况下Nastran的MPC等约束关系相关关键词后面第一个参数EID是Element Number，如下图，所以猜测采用的单元形式。

4.2 Abaqus的实现方式

Abaqus的关键词中inp的MPC中已经消除了这个单元号，如下图，Abaqus如果是按单元来实现的话应该会在inp关键词中反应，譬如质量点，虽然Abaqus/CAE上是inertia设置，但inp依然是Element关键词，也会设置单元类型为MASS，但MPC猜测没有这么做，而是把Master点和Slave点的关系写到了关系矩阵中，在组装全局刚度阵时将关系式单独写到全局刚度阵中。

4.3 iSolver的实现方式

iSolver实现时一开始是按第二种方法，也就是Abaqus的方式来做的，但发现这种方法当遇到复杂问题时改动比较大，无论哪种方法，Master节点自由度和Slave之间的关系都是局部节点的系数，这个系数在全局节点下都需要进一步转换，而第二种方法就需要再自己转换一次，当遇到复杂问题譬如两个约束关系嵌套（一个约束关系的Master是另一个约束关系的Slave）、或者一个Slave被两个Master用到时这个转换的编程复杂度也会上去，后面我们对于部分约束关系改为了单元形式，因为此时就只需要在整体组装时一次性的将局部节点自由度关系转换到全局关系了。

如果你自己编程序，不管采用哪种实现方式，我们觉得除了要实现这个功能外，还要考虑你的程序的兼容性，不破坏你原有的架构。

软件具体实现的问题就不再详述了，我们再接着讨论两个Lagrange因子求解约束方程两个基本的容易搞混的问题：

5 约束方程的个数

约束关系具体是各个节点自由度之间的关系，但为了简单起见，我们只认为是节点之间的约束关系，抛开自由度个数，也就是无论节点是3个还是6个，只要还是这个节点的约束方程，个数只认为是1个。

那么每种约束关系（每个节点下的具体的每个自由度约束对应的一个Lagrange因子）有多少个约束方程？

1） 点之间的约束关系

RBE2每个Slave节点的位移都必须和Master节点相等，也就是每个Slave节点对应1个约束关系，所以RBE2的约束关系数目=Slave节点数，对每一个Slave节点自由度，都应该有一个Lagrange因子，该Lagrange因子可以认为绑定在Slave节点上的物理量。而RBE3中因为所有的Slave和一个Master节点组成了一个约束关系，所以只有一个真正的约束关系，所以RBE3的约束关系数目=Master节点数（一般是1个），Lagrange因子可以认为绑定在Master节点。

2）面之间的约束关系

由面之间的约束关系方程，可知对每个Slave节点，都有一个约束方程，因此， Contact或者Tie等面之间的约束方程个数=Slave节点数目，和RBE2一致。Lagrange因子可以认为绑定在Slave节点。

6 Lagrange因子的物理含义

由上述方程从形式上来看，可知，Lagrange因子类比普通单元中的应力， 而h类比普通单元中的应变，但还是很多人无法理解Lagrange因子到底是什么，我们理解也是实际的物理量。分两种情况

6.1 点之间约束的Lagrange因子含义

对于点之间的约束关系的h为Master点和Slave点的位移差（简单起见，我们不考虑转动），单位是长度

6.2 面之间约束的Lagrange因子含义

对于面之间的约束关系的h为Master面和Slave面的方向距离，单位是长度

约束关系的能量项为：

6.3 点之间约束的Lagrange因子含义证明

为证明点约束的Lagrange因子是集中力，我们取一个体单元43.75X10X5及一个RP点，RP点添加x方向集中力载荷1e9。

RP点和体单元采用KCoupling绑定。

6.3.1 单节点绑定

第一次仅绑定8号节点。

采用iSolver求解Lagrange因子得到如下（此时没有打开几何非线性）：

可以发现8号节点Lagrange因子为-1e9、1.95e-3、1.95e-3，相对e9来说，yz方向的Lagrange因子可以略为0，得到该节点正好和右边的x方向的力平衡，和前面面之间约束关系类似，KCoupling的Lagrange因子为Slave节点处Slave节点对Master节点的作用力，在图上表示如下。方向或者正负号完全由约束方程的正负决定，iSolver是这种负号是因为取的约束方程是h=Us-Um，而不是h=Um-Us。

而且这个作用力在当前时刻猜测不是沿着Slave节点指向Master节点，因为此时iSolver计算出来的节点8的位移在yz方向要明显小于RP点位移，所以，8节点指向RP点已经不是沿x方向了，也就是Slave和Master节点之间不是类似杆只承受拉压，而是类似梁可以承受切向力。

6.3.2 多节点绑定

Slave采用7号和8号节点绑定。如下图所示：

采用iSolver求解Lagrange因子得到：

可以发现7号节点的z正好和8号节点z方向作用力抵消，而8号节点的x方向的作用力正好和外力抵消。

其实KCoupling绑定的所有Slave节点Lagrange因子的和肯定和外力相同，这点可以直接用约束关系推导求出，有兴趣的可以试一下。

6.3.3 问题

上述点之间约束关系的Lagrange是Slave节点和Master节点之间的作用力怎么用Abaqus来验证一直没有找到类似接触的CPRESS合适的对应物理量，要是哪位大神知道也请告知一下。

6.4 面之间约束的Lagrange因子含义证明

   为证明面约束的Lagrange因子是单位面积上的接触力，我们取两个5X1X0.1的长方体，设置两个面紧贴在一起，设置不滑动，下方Master体单元的四个点固支，上方Slave单元的四个点加法向1e10的压力。

6.4.1 理论值

显然，理论上的面压力=1e10*4/(5X1)=8e9。

6.4.2 iSolver解

在iSolver中，我们求取Lagrange因子采用和Slave节点相同的自由度所在全局中的位置。查询单元，可知Slave单元为1号单元，由5,6,8,7,1,2,4,3组成，而Slave节点为8,7,4,3。每个节点三个自由度，即Lagrange因子所在自由度为7-12，和19-24。

采用iSolver求解器求解Lagrange因子，打印出Lagrange的值：

可发现第9和12的Lagrange因子的值为-8e9，和理论一致。由于是负值，所以指向Z-方向，表示Slave节点处Slave节点对Master面的压力。

6.4.3 Abaqus解

将iSolver得到的模型输出到Abaqus进行求解，暂时没找到Abaqus直接的Lagrange因子打印方法，但既然我们认为Lagrange就是面力，可以查看Abaqus接触面力输出。

Abaqus接触面压力为CPRESS，查看CPRESS如下图所示，为8e9，猜测CPRESS是负的Lagrange因子，表示Slave节点处Master对Slave节点的面压力。

同时，我们也可以查看Abaqus的单位面上的切向摩擦力，即CSHEAR1/2，如下所示：

显然，和iSolver的Lagrange因子正好差个负号，猜测也是取了负的Lagrange因子。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删