有限元基础编程——再论罚函数

引言

之前一期讨论过如何使用罚函数来处理边界条件，但是讲的不具有普遍性，即求解时整体刚度矩阵行列已经被划掉了，然后处理“斜支座”的问题，这次将从系统引入罚函数后的势能函数（泛函）开始讲起，然后求解固支边界条件、固支与位移耦合混合的边界条件，两个方面应用罚函数处理边界条件，在求解时保持整体刚度矩阵大小不变。

系统势能建立

给定位移的边界条件

考虑边界条件    ，其中    是沿支座自由度 1 的已知给定位移。用一个具有较大刚度C的弹簧来模拟支座（C一般取刚度矩阵中最大元素的10000倍），此时弹簧的一端移动了    ，因为结构在该点处产生的附加抗力相当小（荷载列阵中忽略该自由度对应的支反力），故沿自由度 1 的位移    近似等于    ，因此弹簧的实际伸长量为    。弹簧的应变能：   

 整体势能函数为：

对势能取最小值    ：

平衡方程展开为：

 操作细节：将大数    加到    的一个对角元上并将    加到    上，此时的    可看作为外荷载列阵，忽略支反力。节点 1 处的支反力等于弹簧在结构上施加的力，弹簧的实际伸长量为    。

支反力为：   

多点约束与已知约束的混合处理

针对特定情况：斜滚子支座、刚性连接等的边界条件形式为：

其中，    为已知常数。系统总势能为： 

 操作细节：

1. 1. 将罚刚度 C 加到对应已知自由度的对角线整体刚度矩阵中，相应的外载荷列阵做出上节描述中的更新；

2. 2. 开始继续修正整体刚度矩阵、外荷载列阵:

3.      
   

3. 支反力计算

以上概念性的东西讲的7788了，接下来我们通过一个算例来理解上述过程的概念。

算例分析

案例1

为了排版简洁，我们以一个简单案例（《工程中的有限元方法 第三版》P80）来描述罚函数如何处理已知位移边界条件。

系统原整体刚度矩阵为：

整体载荷列阵为：

注意：这里的载荷列阵没有加入支反力影响哦~上述过程与乘大数法极为相似。现在我们开始加入罚刚度：因为自由度 1 和 3 是固定的，故将罚刚度 C 加到整体矩阵第 1 和第 3 对角元上，取 C 为    。

整体刚度平衡方程变为：

求解得出节点位移    ，然后根据上述公式求得支反力。

小结：

1. 1. 刚度矩阵要知道在哪个对角线上加上罚刚度；

2. 2. 平衡方程中的位移列阵中的元素都视为未知求解量；

3. 3. 方程右端荷载列阵中，忽略支反力因素，如果固支边界条件，则为0，如果已知位移为      ，则需要在荷载列阵对应元素将0替换为      。

案例2

参照《工程中的有限元方法 第三版》P83》。

位移边界条件：

系统原整体刚度矩阵：

接下来加入罚刚度：

1. 1. K(3,3)和K(4,4)对角元素加上C；

2. 2. 考虑第一个约束方程组：

3.      
   

4.     将其加入到：

5.    

6. 3. 考虑第二个约束方程组，与上述同样的步骤。

7. 最终得到整体刚度方程：

8.      
   

9. 对上述矩阵方程进行求解即可得到节点位移，然后根据本文公式求得节点支反力即可。

以上即为木木理解的罚函数处理边界条件的”基本套路“，在实际编程的时候可编制相应程序来计算，主要先是理解理论概念，编程这一块属于后话，用什么语言编，简洁与否等等可以后期多练多遍。