将单元刚度矩阵组装为全局刚度矩阵后,有: 此时的线性方程没有唯一解,
.边界条件分为两类:Forced and Geometric;对于力边界条件可以直接附加到节点力向量
因此我们基本只需要处理Geometric Boundary condition.下面介绍三种方法,将Bcs引入到
、
以位移边界条件为例,指定相关自由度值即:
Method 1 将开头的
其中,已知节点力 向量;
公式2进一步:
这时,
一旦
Method 2 也称划行划列法 .method 1 中需要对
实际计算中,不需要对刚度阵重新排序.算法操作如下:
对所有的约束自由度
Method 3 该方法也称乘大数法 .假设约束自由度为
该方法通用性强,适合大多数的静力学线性问题,但数值精度与大数的取值有关,太小了精度差,太大了容易出现"矩阵奇异"的现象
Method 4(对角元素置1法) 该方法的做法是,对于约束自由度
以6x6的刚度矩阵为例子,
不引入大数,避免了数值稳定性的问题,不会影响矩阵的条件数; 但只适合 ;可能影响系统矩阵的特性,直接替换可能改变矩阵的对称性(尤其在动力学和非线性问题中);不能处理非0的位移加载,只能处理力加载
Example 例题来自《The Finite Element Method in Engineering》的悬臂梁模型(example6.4, page227)
静力平衡方程为:
解为:
solve by method 1 solve by method 2 循环每个位移约束,需要注意高亮处的操作 :
求解:
solve by method 3 Code Realize 四种方法进行Python+Numpy+Scipy编程实现,并与Example的解进行对比.
import numpy as npfrom typing import Dict ,List ,Tuple import scipy as scdef Method1 (K:np.ndarray,BcDict:Dict [int ,float ],LoadDict:Dict [int ,float ] )->Tuple :0 ]1 )),np.zeros((dofNum,1 ))list (BcDict.keys()))-1 list (set (range (dofNum))-set (prescribedDofIndexs.tolist())),dtype=int )for label,Pval in LoadDict.items():1 0 ]+=Pvalfor label,Uval in BcDict.items():1 0 ]+=Uval0 ])].copy()0 ])].copy()0 ])]=U10 ])]=P2return U,Pdef Method2 (K:np.ndarray,BcDict:Dict [int ,float ],LoadDict:Dict [int ,float ] )->Tuple :0 ]1 )),np.zeros((dofNum,1 ))for label,Pval in LoadDict.items():1 0 ]+=Pvalfor label,Uval in BcDict.items():1 for i in range (dofNum):0 ]=P[i,0 ]-K[i,j]*Uvalfor i in range (dofNum):0 0 1 0 ]=Uvalreturn U_,P_def Method3 (K:np.ndarray,BcDict:Dict [int ,float ],LoadDict:Dict [int ,float ] )->Tuple :max (K)*10e6 0 ]1 )),np.zeros((dofNum,1 ))for label,Pval in LoadDict.items():1 0 ]+=Pvalfor label,Uval in BcDict.items():1 0 ]=K[j,j]*Uvalreturn U_,P_def Method4 (K:np.ndarray,BcDict:Dict [int ,float ],LoadDict:Dict [int ,float ] )->Tuple :if np.any (np.array(list (BcDict.values())) != 0 ):raise ValueError('该方法不能处理非0位移加载' )0 ]1 )),np.zeros((dofNum,1 ))for label,Pval in LoadDict.items():1 0 ]+=Pvalfor label, Uval in BcDict.items():1 0.0 0.0 1.0 0 ]=0 return U_,P_if __name__ == '__main__' :12 ,600 ,-12 ,600 ],600 ,40000 ,-600 ,20000 ],12 ,-600 ,12 ,-600 ],600 ,20000 ,-600 ,40000 ]])1 :0 ,2 :0 }3 :-50 ,4 :20 }0 ,0 ,-16.5667 ,-0.2480 ])50.0 ,4980.0 ,-50 ,20 ])print (f"extract U=\n{extract_U} " )print (f"u=\n{u.T} " )print (f"extract_P=\n{extract_P} " )print (f"p=\n{p.T} " )计算结果:
extract_U=总结 列举了四种直接节点位移边界条件的处理办法,并编程实现,求解案例.对比结果发现:相比Method3存在数值误差,其他三个都更加精确.
如果需要处理多点耦合边界条件,则有罚函数法,拉格朗日乘子法等 .
参考资料:
Note Completed at 2025/03/08
为何封面二次元，因为我是一个lsp（滑稽.jpg）
