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

技术分享 11个月前 (03-09) 0 999+

关注

将单元刚度矩阵组装为全局刚度矩阵后,有:

此时的线性方程没有唯一解,([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}_1tag{2} ]

[[K_{21}]overrightarrow{Phi}_1+[K_{22}]overrightarrow{Phi}_2=overrightarrow{P}_2tag{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[0right] \ left[0right] & left[Iright] 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} ]

实际计算中,不需要对刚度阵重新排序.算法操作如下:

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

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

Method 3

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

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

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)

静力平衡方程为:

解为:

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

[P=[50.0,4980.0,-50,20] ]

solve by method 1

solve by method 2

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

求解:

solve by method 3

Code Realize

四种方法进行Python+Numpy+Scipy编程实现,并与Example的解进行对比.

#------------------------------------------------------------------------------- # Name:        BcsProcess # Purpose:     引入边界条件到[K]中,并返回解[U],[P] #               input: #                   K:全局刚度矩阵,(M,M) numpy.array #                   BcDict:位移约束,key (int) = 自由度序号(1-based) , value (float) = 自由度约束值 #                   LoadDict:节点力加载,key (int) = 自由度序号(1-based) , value (float) = 施加的节点力加载或者等效节点力加载 # # Author:      Administrator # # Created:     08-03-2025 # Copyright:   (c) Administrator 2025 # Licence:     <your licence> #------------------------------------------------------------------------------- import numpy as np from typing import Dict,List,Tuple import scipy as sc def Method1(K:np.ndarray,BcDict:Dict[int,float],LoadDict:Dict[int,float])->Tuple:     dofNum=K.shape[0]     # 初始化向量     U,P=np.zeros((dofNum,1)),np.zeros((dofNum,1))     prescribedDofIndexs=np.array(list(BcDict.keys()))-1      #使用集合运算,全部自由度与约束自由度求差, 得到自由位移自由度的     freeDofIndexs=np.array(list(set(range(dofNum))-set(prescribedDofIndexs.tolist())),dtype=int)      # 已知节点力加到P     for label,Pval in LoadDict.items():         ind=label-1         P[ind,0]+=Pval     # 已知节点位移(prescribed dof)     for label,Uval in BcDict.items():         ind=label-1         U[ind,0]+=Uval     U2=U[np.ix_(prescribedDofIndexs,[0])].copy()     # 已知节点力(free dof)     P1=P[np.ix_(freeDofIndexs,[0])].copy()     # 重新划分K行列     K11=K[np.ix_(freeDofIndexs,freeDofIndexs)].copy()     K12=K[np.ix_(freeDofIndexs,prescribedDofIndexs)].copy()     K21=K[np.ix_(prescribedDofIndexs,freeDofIndexs)].copy()     K22=K[np.ix_(prescribedDofIndexs,prescribedDofIndexs)].copy()     # 计算自由节点位移值     U1=np.dot(sc.linalg.inv(K11),P1-K12.dot(U2))     # 计算支反力     P2=np.dot(K21,U1)+np.dot(K22,U2)      # 合并到U,P向量     U[np.ix_(freeDofIndexs,[0])]=U1     P[np.ix_(prescribedDofIndexs,[0])]=P2     return U,P  def Method2(K:np.ndarray,BcDict:Dict[int,float],LoadDict:Dict[int,float])->Tuple:     K_origin=K.copy()     dofNum=K.shape[0]     # 初始化向量     U,P=np.zeros((dofNum,1)),np.zeros((dofNum,1))     # 已知节点力加到 P     for label,Pval in LoadDict.items():         ind=label-1         P[ind,0]+=Pval     # 循环所有的位移约束     for label,Uval in BcDict.items():         j=label-1         #Step1         for i in range(dofNum):             P[i,0]=P[i,0]-K[i,j]*Uval         #Step2         for i in range(dofNum):             K[i,j]=0             K[j,i]=0         K[j,j]=1         #Step3         P[j,0]=Uval      # 求解 K'U=P'     U_=sc.linalg.solve(K,P)     P_=np.dot(K_origin,U_)      return U_,P_  def Method3(K:np.ndarray,BcDict:Dict[int,float],LoadDict:Dict[int,float])->Tuple:     C=np.max(K)*10e6     K_origin=K.copy()     dofNum=K.shape[0]     # 初始化向量     U,P=np.zeros((dofNum,1)),np.zeros((dofNum,1))     # 已知节点力加到 P     for label,Pval in LoadDict.items():         ind=label-1         P[ind,0]+=Pval     # 循环所有位移约束     for label,Uval in BcDict.items():         j=label-1         # Step1         K[j,j]=K[j,j]*C         # Step2         P[j,0]=K[j,j]*Uval      # 求解 K'U=P'     U_=sc.linalg.solve(K,P)     P_=np.dot(K_origin,U_)      return 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位移加载')     K_origin=K.copy()     dofNum=K.shape[0]     # 初始化向量     U,P=np.zeros((dofNum,1)),np.zeros((dofNum,1))      # 已知节点力加到 P     for label,Pval in LoadDict.items():         ind=label-1         P[ind,0]+=Pval      # loop all nodal bcs     for label, Uval in BcDict.items():         j=label-1         K[j,:]=0.0         K[:,j]=0.0         K[j,j]=1.0         P[j,0]=0      # solve K'U=P'     U_=sc.linalg.solve(K,P)     P_=np.dot(K_origin,U_)      return U_,P_  if __name__ == '__main__':     K=np.array([[12,600,-12,600],                 [600,40000,-600,20000],                 [-12,-600,12,-600],                 [600,20000,-600,40000]])     Bcs={1:0,2:0}     loads={3:-50,4:20}      # 精确解     extract_U=np.array([0,0,-16.5667,-0.2480])     extract_P=np.array([50.0,4980.0,-50,20])     # 求解     u,p=Method4(K,Bcs,loads)      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= [  0.       0.     -16.5667  -0.248 ] extract_P= [  50. 4980.  -50.   20.]  solving by method 1 u= [[  0.           0.         -16.56666667  -0.248     ]] p= [[  50. 4980.  -50.   20.]]  solving by method 2 u= [[  0.           0.         -16.56666667  -0.248     ]] p= [[  50. 4980.  -50.   20.]]  solving by method 3 u= [[-1.04166667e-11 -3.11250000e-13 -1.65666667e+01 -2.48000000e-01]] p= [[  50. 4980.  -50.   20.]]  solving by method 4 u= [[  0.           0.         -16.56666667  -0.248     ]] p= [[  50. 4980.  -50.   20.]]