有限元约束关系 | Lagrange因子法求解

3月前浏览4016

本次主要分享内容：

概述
约束关系
Lagrange因子法求解约束方程
各个软件的约束关系的具体实现方式
4.1 Nastran的实现方式
4.2 Abaqus的实现方式
4.3 iSolver的实现方式
约束方程的个数
Lagrange因子的物理含义
6.1 点之间的约束关系
6.2 面之间的约束关系
6.3 点之间约束的Lagrange因子含义证明
6.3.1 单节点绑定
6.3.2 多节点绑定
6.3.3 问题
6.4 面之间约束的Lagrange因子含义证明
6.4.1 理论值
6.4.2 iSolver解
6.4.3 Abaqus解

概述

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

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

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

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

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

约束关系

我们在系列文章第35章介绍了接触分析及常用的基本算法。现在Abaqus、LS-DYNA、Ansys等结构商软都说可以处理复杂的上万零部件接触的整车、整机等模型仿真，没做过实际的这种仿真分析，很好奇，接触分析算法往往涉及大变形、边界不连续，只要输入条件或者算法稍微变化一些，两个零部件算出来的接触结果就可能差异很大，更不用说上万个零部件的接触结果了，对这种大规模组装模型的仿真结果不知如何来判断它的可靠性，像普通的只校核一下材料的应力还是看一下动画是否和试验一致？毕竟仿真只有简单的标准来判断结果的正确性才能在企业中起到真正辅助设计的作用，如果你恰好做过，不妨也简单介绍一下你的经验。对自研CAE软件开发者来说，自研结构CAE软件是否真的要和商软去比拼接触等复杂算法还是花更多时间在精雕细琢那些常用功能上，这也是开发者需要慎重考虑的问题，而且很多自主CAE软件连常规线性问题都算的不对，或者都没法用鼠标稳定的走完那些材料、属性、边界、加载等流程，用户又怎么会相信你能算对接触这种复杂问题的？

不管怎么样，从有限元实现的角度来讲，如果想做真正实际工程中的接触分析，那么首先需要去研究约束关系，接触分析在有限元中也仅是约束关系的一种。有限元中的约束很多场景大家用的是边界中的简支、固支等约束，但从更广泛的角度上讲，只要表示一个节点的某个自由度依赖于其它的节点自由度或者取某个特定值，就可以称为约束关系。只不过对固支、简支等直接自由度=0，在有限元中直接减缩刚度阵就行，很容易求，但对节点自由度相互依赖的约束关系就比较复杂了。约束关系主要有两类。

一类是MPC点之间的约束。Nastran的MPC的灵活度要远远超过Abaqus，Nastran的主节点可以选择123自由度，也可以对每个从节点设置不同的自由度，还能主节点和从节点互相包含，Abaqus更多的是只负责80%的常用应用场景，复杂功能让你编子程序，但事实上一线仿真工程师又有多少人愿意编子程序呢？这种做法导致虽然Abaqus无论从用户体验、非线性还是商业化都比Nastran好很多，但很多线性的工程复杂问题还是没法替代Nastran。
另一类是Contact、Tie等的面之间的约束关系。在这方面Abaqus要明显强于Nastran了。

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

Lagrange因子法求解约束方程

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

由上一章，去掉所有的积分形式，约束方程可以统一写成如下形式:

对也去掉积分形式，且将体力和面力统一写成，加入约束方程后总的虚功原理统一表示为

因为虚位移是任意的，所以只能虚位移前面的系数=0，得到方程组：

对于一般的Newton迭代法，如结构平衡方程写为

对方程设已知近似根，将函数在点采用泰勒展开，有

即：

上述约束方程组带入后得到

其中：

上式中外力随位移的变化一般不考虑，Abaqus默认也如此，如需考虑，需要在Abaqus加载界面上直接勾选Follow Node Rotation等。
一般都是0。

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

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

既然上一章说了约束关系可以类比普通的有限单元，那么可以把约束关系按单元来实现。
约束关系也就是Master节点自由度和Slave节点自由度之间的系数，直接把这个约束关系写下来放到全局刚度阵就行。

Nastran的实现方式

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

Abaqus的实现方式

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

iSolver的实现方式

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

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

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

约束方程的个数

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

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

点之间的约束关系

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

面之间的约束关系

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

Lagrange因子的物理含义

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

点之间约束的Lagrange因子含义

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

乘以得到能量项：

从左右单位要相同的角度，因为上面没有积分，所以必然是集中力。

面之间约束的Lagrange因子含义

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

约束关系的能量项为：

是接触面的法向，且在面内积分，所以Lagrange因子是面平均的位移法向的力，那么猜测就是压力。

如果存在相对面的切向滑动，那么需要考虑存在摩擦力，约束关系的能量项为：

和分别表示两个切向的相对位移，则同理，和分别表示两个切向的单位面积的摩擦力。注意，s1和s2不一定是两个面相对滑动的位移，不滑动时，摩擦力就是静摩擦力。

而

表示单位面积上的三维正交方向的接触力。单位面积上的力不就是应力吗？这也是前面论述的Lagrange因子可以类比普通单元的应力的原因。

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

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

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

单节点绑定

第一次仅绑定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。