数值模拟基础知识：啥是状态方程

上次说到了罚函数算法 接触基础知识：接触计算中的penalty（罚函数）方法而拉格朗日乘子法（LM）是一种通过引入额外自由度直接满足接触约束的数值算法。其核心思想是将接触力视为独立变量，通过约束强制接触界面非穿透，而非依赖弹簧类比。以下从基本原理、算法实现等方面来说说：一、基本思想与数学模型与penalty将接触力视为一个跟穿透距离有关的弹簧力的思想不同，拉格朗日乘子法通过引入拉格朗日乘子 表示接触面上的法向接触力，将接触条件转化为约束方程： 具体物理解释就是：第一个式子表示两个物体不能互相穿透第二个式子表示当接触的时候 , 而在没有接触的时候.我们可以非常简单的与罚函数（penalty）方法建立联系，那就是： 二、数值实现特征 一般来说对于点对面的方法有两种实现方式，一种是在单元节点上进行检测、另外一种是在积分点处进行检测，就以最简单的节点处进行检测为例，对于从面的节点，进行该节点到主面的投影（住从面的网格大概率不一致），得到投影点的坐标，进而反向求得该点的参数坐标，那么穿透量就可以由以下的公式计算得到 由于一半情况下拉格朗日算法都是在隐式分析中使用，所以要推导其切线刚度矩阵与内力向量，具体方法为，对于上述的约束方程取变分（不懂变分，可以就先理解为微分，然后看[[基础知识：（偏导、方向导数）（变分、微分、线性化）]]） 由于这里我们的 与 被视为两套独立的物理量，其中 针对点对面接触，每一个点有一个“接触力”，所以每一个接触从点就对应一个 ，而 的变分正是上面节点位置的线性组合的变分： 这样我们把这些变分对应的矩阵加到总体矩阵中就会出现如下的形式： 其中 的计算无需赘述，和这个拉格朗日算子法计算接触没关系，就正常的单元计算，而 也就对应着上述 的表达形式，这里 这样如此这般通过上述引入的新的自由度 与引入的约束方程我们就建立了接触力与主面从面相应节点的位移之间的关系，不如penalty方法那么直接，但是也很好理解。对于内力向量也非常直接 物理上可以理解为 是接触力，通过形函数以一定的权重分布施加给这些主从面的节点上。也很合理。上面说的是法向的接触力，其实切向一样，可以引入另外一个变量代表切向的接触力， 不过不同的点在于这里切向的力未必是对的，只是先“试一下”，所以叫做试切向力，因为其并不一定满足摩擦力的基本原则： 其中 为摩擦因子。这里补充说一句，摩擦老有意思了，如果展开讲我觉得能写一大堆。 这里可以看到拉格朗日算子这里没有近似，得出的结果就是完美满足接触的约束条件的，那么为啥很多人不愿意用这个方法呢？有些人说算法写起来复杂，但其实上边的原理推导起来并不复杂，甚至在很多情况下推导还更方便，但实现起来主要有以下几个难点：由于接触需要动态的检测接触位置，每一个迭代接触发生的位置都可能不同，所以 的数量可能就不同，对于有限元矩阵组装提了一些要求大家看上边那个矩阵形式，由于对角线现在出现了一大堆0，原本正定的矩阵现在已经变得不正定了，你的矩阵求解器如果只支持正定矩阵求解，那就解不出来额外的自由度多的话，求解可能慢一些，不过现今的技术，这个不太成问题（老的教科书会特别强调这点，别信！）最后总结一下：LM算法对比罚函数法优势 局限性 无穿透 ：严格满足几何约束 自由度数量补丁 高精度 ：接触力直接与乘子关联 矩阵变为非正定 稳定性强 ：尤其适合软硬结合的接触，不太需要调参 求解规模扩大 既然说了这两个，其实在实际工程中，常结合增广拉格朗日法（Augmented Lagrangian Method）可以认为是penalty与LM方法的结合，请听下回分解。来源：大狗子说数值模拟