有限元分析-破解起重机未知节点位移、应变及支座反力谜团

import numpy as npimport pandas as pdimport torchimport mathimport matplotlib.pyplot as plt%matplotlib inlinenod=pd.read_excel('./input1.xlsx', sheet_name='Sheet1', header=None) #导入input数据，每列所代表的含义，前两列为节点坐标，之后X方向的约束，Y方向的约束，X方向的力，Y方向的力material=pd.read_excel('./input1.xlsx', sheet_name='Sheet2', header=None)#前两列单元左右节点编号，杨氏模量，泊松比#单元和节点都导入完毕materialmaterial = np.array(material)nod = np.array(nod)#将单元和节点函数转换为array格式print(material)#单元刚度矩阵已经建立完毕#装配成总刚n=23#节点数ne=43#单元数K=np.zeros([2*n,2*n])element_stiffness = np.zeros([ne, 4, 4])for m in range(ne):    i=int(material[m,0])#第i个单元左、右节点编号j=int(material[m,1])    xi=nod[i-1,0]    yi=nod[i-1,1]    xj=nod[j-1,0]    yj=nod[j-1,1]#第i个单元左右节点的x、y坐标L=math.sqrt(pow((xj-xi),2)+pow((yj-yi),2))#计算第m个杆单元长度# print(L)s=(yj-yi)/L#杆单元sin值c=(xj-xi)/L#杆单元cos值ea=material[m,2]*material[m,3]/L#单元刚度矩阵前系数AE/Lke=ea*torch.tensor([[c*c,c*s,-c*c,-c*s],[s*c,s*s,-s*c,-s*s],[-c*c,-c*s,c*c,c*s],[-c*s,-s*s,s*c,s*s]])#计算每个单元刚度矩阵ke = np.array(ke)    element_stiffness[m,:,:]=ke #一个杆单元在整体坐标下的单元刚度矩阵，下面四行代码是将所有杆变成在整体坐标下的单元刚度矩阵 :代表所有行和所有列K[2*i-1-1:2*i,2*i-1-1:2*i]=K[2*i-1-1:2*i,2*i-1-1:2*i]+ke[0:2,0:2]    K[2*i-1-1:2*i,2*j-1-1:2*j]=K[2*i-1-1:2*i,2*j-1-1:2*j]+ke[0:2,2:4]    K[2*j-1-1:2*j,2*i-1-1:2*i]=K[2*j-1-1:2*j,2*i-1-1:2*i]+ke[2:4,0:2]    K[2*j-1-1:2*j,2*j-1-1:2*j]=K[2*j-1-1:2*j,2*j-1-1:2*j]+ke[2:4,2:4]Stiffness_save=np.zeros((2*n,2*n))Stiffness_save[:]=K#储存全局刚度矩阵,K的变化不能令其更改Stiffness_save的值Stiffness_saveelement_stiffnessKF = np.zeros([2*n, 1]).astype(np.float64) #边界条件处理，采用对角元乘大数法for lop in range(n):    F[2*lop,0]=F[2*lop,0]+nod[lop,4]    F[2*lop+1,0]=F[2*lop+1,0]+nod[lop,5]    if nod[lop,2]>=1:#如果该节点x方向有约束       K[2*lop,2*lop]=100000000*K[2*lop,2*lop]    if nod[lop,3]>=1:#如果该节点y方向有约束       K[2*lop+1,2*lop+1]=100000000*K[2*lop+1,2*lop+1]Fk_inverse =np.linalg.inv(K)#调用numpy进行函数的矩阵的求逆# print(k_inverse)u=np.matmul(k_inverse,F)#矩阵的乘法Force=np.matmul(Stiffness_save,u)#求出节点支反力force_of_support_nodes=Force[0:4]force_of_support_nodes#查看节点支反力k_inverse[0,0]u=np.matmul(k_inverse,F)#求出节点位移向量Force=np.matmul(Stiffness_save,u)#求出节点支反力force_of_support_nodes=Force[0:4]force_of_support_nodes#查看节点支反力uelement_force=np.zeros((ne,1))#储存单元内力element_stress=np.zeros((ne,1))#储存单元应力for lop in range(n):    i=int(material[lop,0])    j=int(material[lop,1])    E=material[lop,2]#第lop个杆单元杨氏模量A=material[lop,3]#第lop个杆单元截面积xi=nod[i-1,0]    yi=nod[i-1,1]    xj=nod[j-1,0]    yj=nod[j-1,1]    L=math.sqrt(pow((xj-xi),2)+pow((yj-yi),2))#计算第lop个杆单元长度s=(yj-yi)/L#杆单元sin值c=(xj-xi)/L#杆单元cos值element_force[lop,0]=E*A/L*((u[2*j-2,0]-u[2*i-2,0])*c+(u[2*j-1,0]-u[2*i-1,0])*s)#求杆单元内力element_stress[lop,0]=E*((u[2*j-2,0]-u[2*i-2,0])*c+(u[2*j-1,0]-u[2*i-1,0])*s)/L#求杆单元应力max_tensile_stress=np.amax(element_stress)max_compressive_stress=np.amax(element_stress)max_tensile_stress,max_compressive_stress#输出最大拉、压应力u_yaxis=np.zeros(n)#储存竖向位移for i in range(n):    u_yaxis[i]=u[2*i+1]max_u_vertical=np.min(u_yaxis)#得到最大竖向位移row=np.argmin(u_yaxis)#返回最大竖向位移位置max_u_X=nod[row,0]max_u_Y=nod[row,1]max_u_vertical,max_u_X,max_u_Y#查看最大竖向位移节点的x,y坐标#将上述结果写入csv文档r1=[max_u_vertical,max_u_X,max_u_Y]r2=[max_tensile_stress,max_compressive_stress]r3=force_of_support_nodesprint(r1,r2,r3)import csvwith open("Output.csv", "w") as csvfile:    writer = csv.writer(csvfile,lineterminator='\n')    writer.writerow(["最大竖向位移","最大竖向位移x坐标","最大竖向位移y坐标"])    writer.writerows([r1])    writer.writerow(["最大拉应力","最大压应力"])    writer.writerows([r2])    writer.writerow(["1节点x方向反力","1节点y方向反力","2节点x方向反力","2节点y方向反力"])    writer.writerows([r3])max_u_vertical#绘制变形前桁架X=np.zeros(2)Y=np.zeros(2)X_after=np.zeros(2)Y_after=np.zeros(2)plt.figure(figsize=(6,9))print(type(material))for lop in range(ne):    i=int(material[lop,0])#第lop个单元左、右节点编号    j=int(material[lop,1])    X[0],X[1]=nod[i-1,0],nod[j-1,0]    Y[0],Y[1]=nod[i-1,1],nod[j-1,1]    X_after[0]=X[0]+u[2*i-2,0]    X_after[1]=X[1]+u[2*j-2,0]    Y_after[0]=Y[0]+u[2*i-1,0]    Y_after[1]=Y[1]+u[2*j-1,0]    plt.plot(X,Y,color='b',linewidth=2.5,marker='o',markerfacecolor='w',markeredgecolor='b')    plt.plot(X_after,Y_after,color='g',linewidth=2.5,linestyle='--',marker='o',markerfacecolor='r',markeredgecolor='y')plt.savefig("./Output_deformed_shape.jpg")plt.show()