有限元基础编程 | 罚函数处理边界条件（理论推导+数值实现）

系统势能建立

给定位移的边界条件

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

整体势能函数为：

对势能取最小值：

平衡方程展开为：

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

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

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

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

操作细节：

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

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

• 支反力计算

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

算例分析

案例1

《工程中的有限元方法 第三版》P80

系统原整体刚度矩阵为：

整体载荷列阵为：

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

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

小结：

• 刚度矩阵要知道在哪个对角线上加上罚刚度；
• 平衡方程中的位移列阵中的元素都视为未知求解量；
• 方程右端荷载列阵中，忽略支反力因素，如果固支边界条件，则为0，如果已知位移为，则需要在荷载列阵对应元素将0替换为。

案例2

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

位移边界条件：

系统原整体刚度矩阵：

接下来加入罚刚度：

• K(3,3)和K(4,4)对角元素加上C；
• 考虑第一个约束方程组：

  将其加入到：
• 考虑第二个约束方程组，与上述同样的步骤。最终得到整体刚度方程：

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

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

为了方便读者使用，木木给出两版代码，分别用matlab和python进行施加罚刚度，摘自木木自研的MFEA程序片段。

Example Code

使用罚函数处理边界条件时，在求出节点位移后，顺道可以把节点反力求出，非常方便。

Python

# Apply constraints (penalty Method)
penalty = np.max(np.abs(KKG)) * 1e4
for i in range(nconstrain):
    # Check for valid constraint degrees of freedom
    if constrain[i, 1] < 0 or constrain[i, 1] > 2:
        raise ValueError("Constraint degrees of freedom must be within the range 1-3 (adjusted to 0-2 after indexing).")
    m = int((constrain[i, 0] - 1) * Dof + int(constrain[i, 1]))
    KKG[m, m] = KKG[m, m] + penalty
    FFG[m] = FFG[m] + penalty * constrain[i, 2]

# Solve for the nodal displacements
KKG = csr_matrix(KKG)  # Convert the stiffness matrix to CSR format for efficient solving
UUG = spsolve(KKG, FFG)  # Solve the linear system KKG * UUG = FFG

# Solve for the nodal RF
for i in range(nconstrain):
    m = int((constrain[i, 0] - 1) * Dof + int(constrain[i, 1]))
    RFG[m] = penalty * (constrain[i, 2] - UUG[m])

Matlab

penalty = max(abs(KKG(:))) * 1e4;

for i = 1:nconstrain
   m = (obj.Constrain(i,1) - 1) * Dof + obj.Constrain(i,2);
   KKG(m, m) = KKG(m, m) + penalty;
   FFG(m) = FFG(m) + penalty * obj.Constrain(i,3);
end

U = KKG \ FFG;

for i = 1:nconstrain
   m = (obj.Constrain(i,1) - 1) * Dof + obj.Constrain(i,2);
   RF(m) = penalty * (obj.Constrain(i,3) - U(m));
end

现在MFEA程序的三维部分已全面转为Python，当然是按照初版Matlab MFEA的风格来改写的，加入星球（后台发送：星球）的小伙伴可以下载到两版代码，分别是Python和Matlab。

Python版MFEA运行截图

Matlab版MFEA截图