有限元基础编程 | 罚函数处理边界条件（理论推导+数值实现）
系统势能建立
给定位移的边界条件
考虑边界条件，其中是沿支座自由度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截图
