有限元编程中关于边界条件的处理

Rainyboy

引言
在有限编程的一般流程中，在获得原始总刚度矩阵之后，就是引入边界条件了。力边界条件直接表示在激振力向量中，而对于位移边界条件的处理较为复杂。位移边界条件的处理方法有万能钥匙型的“罚函数法”和“拉格朗日乘子法”，以及处理特定情形的“置大数法”和“划0置1”法。而且对于静力学分析，模态分析，瞬态响应分析和谐响应分析，在一些细节上的处理稍有不同。在编程中，可能遇到的技术难点是：1）从理论上，应该选择哪种技术来处理边界条件按更合适？2）在技术细节上，边界条件的处理应该在前处理器中实现，还是在求解器中实现？在这里我打算根据个人的一些经验谈一谈这些问题，抛砖引玉。
1， 位移边界条件的一般形式和“万能钥匙型”处理方法
位移边界条件的一般形式可写为：
{\bf A} {\bf C}={\bf B}
其中{\bf A}和{\bf B}是常数矩阵，而{\bf u}是有限元模型节点自由度列向量。实际上，要解决的问题是，如何在原始动力矩阵中引入这些补充方程的影响，即联解方程组：
{\bf M}{\bf \ddot u} + {\bf C}{\bf \dot u} + {\bf K}{\bf u} = {\bf f}\\
{\bf A}{\bf u}={\bf B}
由于多了若干个（设为m个）位移边界条件，因此上述方程组有N个未知数却有N+m个方程，无法直接联立求解。既然称这里讨论的边界条件形式为“一般形式”，那么能处理这种“一般形式”的方法在这里就姑且称为“万能钥匙型”处理方法。
既然问题的关键在于未知数个数和方程数目的不匹配，那么解决方案也就是如下两种：减少方程个数（罚函数法）和增加未知数个数（拉格朗日乘子法）。在这里不打算展开方法的细节，要知道细节google即可。只是想说明两个方法的区别。罚函数法最终将m个约束方程修正到了原始刚度矩阵中，从而保证了自由度数的不变，这可以视作它的一个优点。但它的缺点也是显而易见的，即需要指定惩罚系数，而且在有多个约束方程的场合，甚至有必要为不同的约束方程指定不同的惩罚系数。而惩罚系数的确定完全是经验的，例如在笔者编写的有限元程序在用罚函数方法来处理压电结构的约束时，需要分别指定针对电场自由度和机械场自由度的惩罚系数才可以达到较满意的精度。拉格朗日乘子法增加了m个未知数，不需要经验的惩罚系数，却增大系统的自由度，在约束方程数量比较多的场合并不令人满意。
注：在ANSYS中，对应于CE命令。

2，特殊位移边界条件及其处理方法：直接位移约束
对于形如
u_{i}=b_{i}\
的位移边界条件，即直接给出部分自由度位移的情形，常见的有两种解决方案，都可以从罚函数方法导出，称为“置大数法”和“划0置1法”。相比较而言，前者所谓的“大数”就是罚函数法中的惩罚系数，因而同样有大数的取值问题。而后者相当于是运用了主对角元素是大数，对方程进行了简化，从而不需要取“大数”，其代价是操作过程稍嫌麻烦。但考虑到“划0置1法”几乎可视作精确地考虑了上述位移边界条件，过程的稍微繁琐似乎也是可以接受的。实际上，若在这两种方法的基础上再往后走，还有一种更“暴力”的方法，姑且称之为“删除行列法”。其基本出发非常好理解，既然一些自由度的取值已经给出，那么它们就可以不作为未知数包含到求解中了。它们的影响体现为一个施加到剩余自由度上的力向量（类似于C-B模态综合法中的约束模态)。其本身对应的行，列将从总体矩阵中删除。这个方法看上去很美，尤其在矩阵中删除行列的操作大多数现代编程语言的数值计算库都能提供的情况下，更显得手到擒来。它唯一的不足是增加了后处理的难度，需要设计数据结构储存那些被删去的自由度的取值并将它们融合到后处理中。
注：在ANSYS中，对应于D命令。

3，特殊位移边界条件及其处理方法：自由度耦合
对于形如
u_{i}=u_{j}=u_{l}=...\

的位移边界条件，通常称为自由度耦合，即表示若干个广义自由度实际上可以由其中任意一个表出。作为一般形式边界条件的一个特例，它当然可以用“万能钥匙型”方法来解决，这里还有一个帖子详细进行了详细的讨论请教结构的自由度耦合的问题。
但实际上我发现有一种更简单直接的方法，既然有几个自由度可以由同一个自由度表出，何不就将它们视为一个自由度？具体地，在内部的“自由度”-“方程编号”映射表中，可以作如下处理：
节点  自由度    方程编号
1        UX         10
6        UX         10
19      UX         10
这样自由度1UX,6Ux和19UX就都被u_{10}表出了。这样做的好处是，所有与这些自由度相关的系数（单元刚度系数，质量系数等）都会组集到同一个实际自由度u_{10}上。当然，其代价是要重新调整“自由度”-“方程编号”映射表，甚至要重新组集总体刚度矩阵。其优点也是显然的，即无需任何经验系数，不增加反而减少了系统的自由度。也许文献中对这种处理方法有命名，有知道的同志请指点一下。这个思路类似于上面的“删除行列法”但是实际操作却并是删除那么简单。
注：在ANSYS中，对应于CP命令。

4,  实现的问题：前处理器还是求解器？
那么，在实际编程中，一个即将遭遇的问题是，对约束条件的处理应当在前处理器中实现还是在求解器中实现？也许根据程序流程思路的不同，你的代码并没有明确的区分前处理器和求解器。但是一个显而易见的问题是代码的复用性，如果边界条件的引入过程都是一样的，为什么不在前处理器中进行，这样，求解器收到的，都是一个已经考虑了边界条件的模型，直接按照各自的算法求解就可以了，不是更好么？问题在于，并不是所有的边界条件都适合在前处理器中进行引入，有些边界条件的处理随着求解器的不同而不同，有些边界条件的处理又与求解器无关。这似乎说明有必要将一部分边界条件在前处理器中编码，一部分在求解器中编码，以实现代码的最佳复用性。这也是我们在“另一个的问题：既然有一般方法，为什么仍然有必要实现特殊方法？”小节中给出的理由之一。
例如，在静力学分析中，处理由约束u_{i}=1引起的节点力只需要考虑静刚度矩阵{f K}；然而在谐响应分析中，则需要考虑动力刚度{f H} = -omega^2{f M} + jomega{f C} + {f K}。姑且不论还有另一个更大的问题，即谐响应分析中约束表达式u_{i}=1
存在歧义，即无法区分它意味着静约束u_{i}(t)=1还是简谐约束u_{i}(t)=1cdot e^{jomega t}。前者产生一个静变形，不在谐响应分析考察的范畴内，将被处理为u_{i}(t)=0；而后者将产生一个简谐力，将被完整地处理。由于这样的歧义，求解器必须知道约束的意义，与此同时，前处理器在不知道下游的求解器是什么时也无法确定如何处理该约束。

这里列出了几种常规求解情况下的约束处理需要修正的原始有限元数据：约束类型（以ANSYS命令命名）  
                                 实际修改数据         静力学分析需要       模态分析需要       谐响应分析需要                                                                                                               的修改数据             的修改数据           的修改数据
D 型                         刚度，力载荷          刚度，力载荷            刚度                 静约束：刚度      
                                                                                                               动约束：刚度，质量，阻尼，力载荷         
CE 型 （罚函数）     刚度，力载荷          刚度，力载荷           刚度                  静约束：刚度
                                                                                                               动约束：刚度，质量，阻尼，力载荷
CE 型 （拉格朗日）  刚度，力载荷         刚度，力载荷           刚度                  静约束：刚度
                                                                                                              动约束：刚度，质量，力载荷
CP 型                        刚度，质量，阻尼             刚度             刚度   质量          静约束：刚度
                                                                                                               动约束：刚度，质量，阻尼
可见，D型和CE型约束的处理与求解器相关，因为它们都会产生等效力载荷向量，这个力载荷怎么计算，是“静”的还是“动”的都只能由求解器来解释。而CP型约束本身不产生约束力，所以与求解器无关。

5, 另一个的问题：既然有一般方法，为什么仍然有必要实现特殊方法？

这个问题扩展开来，是：既然我们有一般处理方法（罚函数法，拉格朗日乘子），那么何不选择一种实现，然后在前处理中或者求解器中调用，将边界条件引入，然后求解？为什么要将约束按照各种类型分开处理？

原因1：如第4节所言，处于代码维护和重用性的角度考虑，CP型约束适合在前处理器中进行，而D型和CE型适合在求解器中进行。既然分开处理，就可以用各自更适用更简单的处理方式了。原因2：一般处理方法各有弊端，选择罚函数必须冒着求解精度下降的风险，而选择拉格朗日乘子则会增大模型的自由度规模。而将特殊情况用各自的特殊方法处理，可以获得精确解。只有在约束条件不特殊时再用一般方法进行处理。这样可以更好的控制模型的精度和规模。

6, 一种实现方式

综上所示，根据个人的编程经验，推荐这样的实现方式：

D 型，  求解器实现，    划0置1法
CE型， 求解器实现，     罚函数法
CP型， 前处理器实现， “重新组集法”

即尽可能采用精确的处理方式，实在避不过时（只能用CE型表示时），用近似解法。
--END--

wangbr2003

了解下情况。。。

mayaview

位移边界讲解的不错！有没有计划讲讲其他类型的边界呢？比如速度边界，接触，MPC等等，感觉会更有意思。

Rainyboy

mayaview 发表于 2015-4-29 07:15
位移边界讲解的不错！有没有计划讲讲其他类型的边界呢？比如速度边界，接触，MPC等等，感觉会更有意思。

摩擦做过一些，但速度边界等等不熟悉了，兄台来一贴？

mxlzhenzhu

我在一本有限元书籍中看到过，说罚函数精度没有拉格朗日法精度高，但是后者增加了自由度数，可能是不希望发生的。在通用有限元nastran里面，保留了两种方法，我看到有的是可以人工干预的，但是有些问题是只能限制一种方法。

mayaview

Rainyboy 发表于 2015-5-1 02:04
摩擦做过一些，但速度边界等等不熟悉了，兄台来一贴？

兄台都算专业摩擦的了，你对接触都不熟悉的话，那就更不敢班门弄斧了。

动力学边界真没有了解过，一直不太清楚细节。

zhangyang1214

院长，如果自由度耦合的情况，自由度之间是具有平方关系如何处理。如u1^2=u2^2+u3^2