登录/ 注册
 找回密码
 立即注册
关闭

CSDN热搜

程序园 精品问答 技术交流 资源下载
返回列表 发新帖
首页 › 业界区 › 业界 ›Note_Fem边界条件的处理和numpy实现的四种方法 ...

Note_Fem边界条件的处理和numpy实现的四种方法

訾懵 2025-6-4 21:55:18
将单元刚度矩阵组装为全局刚度矩阵后,有:
1.png

此时的线性方程没有唯一解,\([K]\)是奇异矩阵,这是没有引入边界条件,消除刚体位移的原因.
边界条件分为两类:Forced and Geometric;对于力边界条件可以直接附加到节点力向量\([P]\)中,即\(P_j=P_j^{*}\),\(P_j^{*}\)是给定的节点力值.
因此我们基本只需要处理Geometric Boundary condition.下面介绍三种方法,将Bcs引入到\([K]、[P]\)
以位移边界条件为例,指定相关自由度值即:\(\Phi_j=\Phi_j^{*}\)
Method 1

将开头的\([K][\Phi]=[P]\)划分为:

\[\begin{bmatrix}[K_{11}] & [K_{12}] \\[K_{21}] & [K_{22}]\end{bmatrix}\begin{Bmatrix}\overrightarrow{\Phi}_1 \\\overrightarrow{\Phi}_2\end{Bmatrix}=\begin{Bmatrix}\overrightarrow{P}_1 \\\overrightarrow{P}_2\end{Bmatrix}\tag{1}\]
其中,\(\Phi_1\)是未知的自由节点自由度向量(free dofs);\(\Phi_2\)是已知的约束节点自由度值\(\Phi_j^{*}\)向量(specified nodal dof);\(P_1\)是已知节点力向量;\(P_2\)是未知的支反力向量
公式2进一步:

\[[K_{11}]\overrightarrow{\Phi}_1+[K_{12}]\overrightarrow{\Phi}_2=\overrightarrow{P}_1\tag{2}\]

\[[K_{21}]\overrightarrow{\Phi}_1+[K_{22}]\overrightarrow{\Phi}_2=\overrightarrow{P}_2\tag{3}\]
这时,\([K_{11}]\)是非奇异矩阵.因此自由节点自由度(未知节点位移)可求:

\[\overrightarrow{\Phi}_1=[K_{11}]^{-1}(\overrightarrow{P}_1-[K_{12}]\overrightarrow{\Phi}_2)\tag{4}\]
一旦\(\Phi_1\)求得,则未知支反力\(P_2\)可由公式3求得.
Method 2

也称划行划列法.method 1 中需要对\([K] ,[\Phi],[P]\)进行行列对调,重新排序.当出现非0位移边界时,method 1耗时长且需要记录过程,之后还需要恢复刚度矩阵.因此和method 1等效的处理方法是构建下式:

\[\begin{bmatrix}\left[K_{11}\right] & \left[0\right] \\\left[0\right] & \left[I\right]\end{bmatrix}\begin{bmatrix}\overrightarrow{\Phi}_1 \\\overrightarrow{\Phi}_2\end{bmatrix}=\begin{bmatrix}\overrightarrow{P}_1-\left[K_{12}\right]\overrightarrow{\Phi}_2\\{\overrightarrow{\Phi}_2}\end{bmatrix}\tag{5}\]
实际计算中,不需要对刚度阵重新排序.算法操作如下:
2.png

3.png

对所有的约束自由度\(\Phi_j\)重复Step 1~3即可,这种操作能够保持刚度和方程的对称性.
Method 3

该方法也称乘大数法.假设约束自由度为\(\Phi_j=\Phi_j^*\),操作如下:
4.png

该方法通用性强,适合大多数的静力学线性问题,但数值精度与大数的取值有关,太小了精度差,太大了容易出现"矩阵奇异"的现象
Method 4(对角元素置1法)

该方法的做法是,对于约束自由度\(\Phi_j=0\),把\([K]\)的j行j列置0,但对角元素Kjj=1,\([P]\)中对应元素置0.
以6x6的刚度矩阵为例子,

\[\begin{gathered}\begin{bmatrix}k_{11} & k_{12} & 0 & k_{14} & k_{15} & k_{16} \\k_{21} & k_{22} & 0 & k_{24} & k_{25} & k_{26} \\0 & 0 & 1 & 0 & 0 & 0 \\k_{41} & k_{42} & 0 & k_{44} & k_{45} & k_{46} \\k_{51} & k_{52} & 0 & k_{54} & k_{55} & k_{56} \\k_{e1} & k_{e3} & 0 & k_{eA} & k_{e5} & k_{e6}\end{bmatrix}\begin{bmatrix}\delta_1 \\\delta_2 \\\delta_3 \\\delta_4 \\\delta_5 \\\delta_6\end{bmatrix}=\begin{bmatrix}f_1 \\f_2 \\0 \\f_4 \\f_5 \\f_6\end{bmatrix}\end{gathered}\]
不引入大数,避免了数值稳定性的问题,不会影响矩阵的条件数; 但只适合\(\Phi_j=0\)这样的简单边界;可能影响系统矩阵的特性,直接替换可能改变矩阵的对称性(尤其在动力学和非线性问题中);不能处理非0的位移加载,只能处理力加载
Example

例题来自《The Finite Element Method in Engineering》的悬臂梁模型(example6.4, page227)
5.png

静力平衡方程为:
6.png

解为:

\[W=[0,0,-16.5667,-0.2480]\]

\[P=[50.0,4980.0,-50,20]\]
solve by method 1

7.png

8.png

9.png

solve by method 2

10.png

循环每个位移约束,需要注意高亮处的操作:
11.png

求解:
12.png

solve by method 3

13.png

14.png

Code Realize

四种方法进行Python+Numpy+Scipy编程实现,并与Example的解进行对比.
  1. #-------------------------------------------------------------------------------
  2. # Name:        BcsProcess
  3. # Purpose:     引入边界条件到[K]中,并返回解[U],[P]
  4. #               input:
  5. #                   K:全局刚度矩阵,(M,M) numpy.array
  6. #                   BcDict:位移约束,key (int) = 自由度序号(1-based) , value (float) = 自由度约束值
  7. #                   LoadDict:节点力加载,key (int) = 自由度序号(1-based) , value (float) = 施加的节点力加载或者等效节点力加载
  8. #
  9. # Author:      Administrator
  10. #
  11. # Created:     08-03-2025
  12. # Copyright:   (c) Administrator 2025
  13. # Licence:     <your licence>
  14. #-------------------------------------------------------------------------------
  15. import numpy as np
  16. from typing import Dict,List,Tuple
  17. import scipy as sc
  18. def Method1(K:np.ndarray,BcDict:Dict[int,float],LoadDict:Dict[int,float])->Tuple:
  19.     dofNum=K.shape[0]
  20.     # 初始化向量
  21.     U,P=np.zeros((dofNum,1)),np.zeros((dofNum,1))
  22.     prescribedDofIndexs=np.array(list(BcDict.keys()))-1
  24.     #使用集合运算,全部自由度与约束自由度求差, 得到自由位移自由度的
  25.     freeDofIndexs=np.array(list(set(range(dofNum))-set(prescribedDofIndexs.tolist())),dtype=int)
  27.     # 已知节点力加到P
  28.     for label,Pval in LoadDict.items():
  29.         ind=label-1
  30.         P[ind,0]+=Pval
  31.     # 已知节点位移(prescribed dof)
  32.     for label,Uval in BcDict.items():
  33.         ind=label-1
  34.         U[ind,0]+=Uval
  35.     U2=U[np.ix_(prescribedDofIndexs,[0])].copy()
  36.     # 已知节点力(free dof)
  37.     P1=P[np.ix_(freeDofIndexs,[0])].copy()
  38.     # 重新划分K行列
  39.     K11=K[np.ix_(freeDofIndexs,freeDofIndexs)].copy()
  40.     K12=K[np.ix_(freeDofIndexs,prescribedDofIndexs)].copy()
  41.     K21=K[np.ix_(prescribedDofIndexs,freeDofIndexs)].copy()
  42.     K22=K[np.ix_(prescribedDofIndexs,prescribedDofIndexs)].copy()
  43.     # 计算自由节点位移值
  44.     U1=np.dot(sc.linalg.inv(K11),P1-K12.dot(U2))
  45.     # 计算支反力
  46.     P2=np.dot(K21,U1)+np.dot(K22,U2)
  48.     # 合并到U,P向量
  49.     U[np.ix_(freeDofIndexs,[0])]=U1
  50.     P[np.ix_(prescribedDofIndexs,[0])]=P2
  51.     return U,P
  53. def Method2(K:np.ndarray,BcDict:Dict[int,float],LoadDict:Dict[int,float])->Tuple:
  54.     K_origin=K.copy()
  55.     dofNum=K.shape[0]
  56.     # 初始化向量
  57.     U,P=np.zeros((dofNum,1)),np.zeros((dofNum,1))
  58.     # 已知节点力加到 P
  59.     for label,Pval in LoadDict.items():
  60.         ind=label-1
  61.         P[ind,0]+=Pval
  62.     # 循环所有的位移约束
  63.     for label,Uval in BcDict.items():
  64.         j=label-1
  65.         #Step1
  66.         for i in range(dofNum):
  67.             P[i,0]=P[i,0]-K[i,j]*Uval
  68.         #Step2
  69.         for i in range(dofNum):
  70.             K[i,j]=0
  71.             K[j,i]=0
  72.         K[j,j]=1
  73.         #Step3
  74.         P[j,0]=Uval
  76.     # 求解 K'U=P'
  77.     U_=sc.linalg.solve(K,P)
  78.     P_=np.dot(K_origin,U_)
  80.     return U_,P_
  82. def Method3(K:np.ndarray,BcDict:Dict[int,float],LoadDict:Dict[int,float])->Tuple:
  83.     C=np.max(K)*10e6
  84.     K_origin=K.copy()
  85.     dofNum=K.shape[0]
  86.     # 初始化向量
  87.     U,P=np.zeros((dofNum,1)),np.zeros((dofNum,1))
  88.     # 已知节点力加到 P
  89.     for label,Pval in LoadDict.items():
  90.         ind=label-1
  91.         P[ind,0]+=Pval
  92.     # 循环所有位移约束
  93.     for label,Uval in BcDict.items():
  94.         j=label-1
  95.         # Step1
  96.         K[j,j]=K[j,j]*C
  97.         # Step2
  98.         P[j,0]=K[j,j]*Uval
  100.     # 求解 K'U=P'
  101.     U_=sc.linalg.solve(K,P)
  102.     P_=np.dot(K_origin,U_)
  104.     return U_,P_
  106. def Method4(K:np.ndarray,BcDict:Dict[int,float],LoadDict:Dict[int,float])->Tuple:
  107.     if np.any(np.array(list(BcDict.values())) != 0):
  108.        raise ValueError('该方法不能处理非0位移加载')
  109.     K_origin=K.copy()
  110.     dofNum=K.shape[0]
  111.     # 初始化向量
  112.     U,P=np.zeros((dofNum,1)),np.zeros((dofNum,1))
  114.     # 已知节点力加到 P
  115.     for label,Pval in LoadDict.items():
  116.         ind=label-1
  117.         P[ind,0]+=Pval
  119.     # loop all nodal bcs
  120.     for label, Uval in BcDict.items():
  121.         j=label-1
  122.         K[j,:]=0.0
  123.         K[:,j]=0.0
  124.         K[j,j]=1.0
  125.         P[j,0]=0
  127.     # solve K'U=P'
  128.     U_=sc.linalg.solve(K,P)
  129.     P_=np.dot(K_origin,U_)
  131.     return U_,P_
  133. if __name__ == '__main__':
  134.     K=np.array([[12,600,-12,600],
  135.                 [600,40000,-600,20000],
  136.                 [-12,-600,12,-600],
  137.                 [600,20000,-600,40000]])
  138.     Bcs={1:0,2:0}
  139.     loads={3:-50,4:20}
  141.     # 精确解
  142.     extract_U=np.array([0,0,-16.5667,-0.2480])
  143.     extract_P=np.array([50.0,4980.0,-50,20])
  144.     # 求解
  145.     u,p=Method4(K,Bcs,loads)
  147.     print(f"extract U=\n{extract_U}")
  148.     print(f"u=\n{u.T}")
  149.     print(f"extract_P=\n{extract_P}")
  150.     print(f"p=\n{p.T}")
复制代码
计算结果:
  1. extract_U=
  2. [  0.       0.     -16.5667  -0.248 ]
  3. extract_P=
  4. [  50. 4980.  -50.   20.]
  6. solving by method 1
  7. u=
  8. [[  0.           0.         -16.56666667  -0.248     ]]
  9. p=
  10. [[  50. 4980.  -50.   20.]]
  12. solving by method 2
  13. u=
  14. [[  0.           0.         -16.56666667  -0.248     ]]
  15. p=
  16. [[  50. 4980.  -50.   20.]]
  18. solving by method 3
  19. u=
  20. [[-1.04166667e-11 -3.11250000e-13 -1.65666667e+01 -2.48000000e-01]]
  21. p=
  22. [[  50. 4980.  -50.   20.]]
  24. solving by method 4
  25. u=
  26. [[  0.           0.         -16.56666667  -0.248     ]]
  27. p=
  28. [[  50. 4980.  -50.   20.]]
复制代码
总结

列举了四种直接节点位移边界条件的处理办法,并编程实现,求解案例.对比结果发现:相比Method3存在数值误差,其他三个都更加精确.
如果需要处理多点耦合边界条件,则有罚函数法,拉格朗日乘子法等.
参考资料:

  • 有限元基础编程 | 边界条件专题(对角元素置"1"法、乘大数法、划行划列法、拉格朗日乘子法、罚函数法)
  • 范雨有限元博客
  • 有限元软件开发  致力于国产大型通用商业有限元计算软件的开发
  • Singiresu S. Rao, sixth edition
Note Completed at 2025/03/08

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！

使用道具 举报
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册
发帖
訾懵
2025-6-4 21:55:18

0

粉丝关注

20

主题发布

板块介绍填写区域，请于后台编辑
微信扫一扫