用python进行简单有限元分析——二维4单元9节点总刚度组装以及9个节点的力和位移求解

2020年数值分析作业，已成功实现，可直接复制代码运行！！

1、大概的理论。这里主要讲总刚度矩阵的组装原理，如何得到单元刚度矩阵请看添加链接描述https://blog.csdn.net/Youngist/article/details/106651143

2、实现得到4单元9节点总刚度矩阵的源代码

from math import *
import numpy as np
from matplotlib import pyplot as plt# 1个单元四个节点
def shapefunction(r,s):#形函数N1 = 1 / 4 * (1 - r) * (1 - s)N2 = 1 / 4 * (1 + r) * (1 - s)N3 = 1 / 4 * (1 + r) * (1 + s)N4 = 1 / 4 * (1 - r) * (1 + s)return N1,N2,N3,N4def diffNdr(r,s): # 求dNidrdN1dr = 1 / 4 * (-1) * (1 - s)dN2dr = 1 / 4 * (1) * (1 - s)dN3dr = 1 / 4 * (1) * (1 + s)dN4dr = 1 / 4 * (-1) * (1 + s)dNdr = [dN1dr,dN2dr,dN3dr,dN4dr]return dNdrdef diffNds(r,s): # 求dNidsdN1ds = 1 / 4 * (1 - r) * (-1)dN2ds = 1 / 4 * (1 + r) * (-1)dN3ds = 1 / 4 * (1 + r) * (1)dN4ds = 1 / 4 * (1 - r) * (1)dNds = [dN1ds, dN2ds, dN3ds, dN4ds]return dNdsdef jacobian(x,y,r,s): # 求J,Jinv,JdetdNdr = diffNdr(r,s)dNds = diffNds(r,s)dxdr = x[0]*dNdr[0]+x[1]*dNdr[1]+x[2]*dNdr[2]+x[3]*dNdr[3]dxds = x[0]*dNds[0]+x[1]*dNds[1]+x[2]*dNds[2]+x[3]*dNds[3]dydr = y[0]*dNdr[0]+y[1]*dNdr[1]+y[2]*dNdr[2]+y[3]*dNdr[3]dyds = y[0]*dNds[0]+y[1]*dNds[1]+y[2]*dNds[2]+y[3]*dNds[3]J = np.array([[dxdr,dxds],[dydr,dyds]])Jdet = np.linalg.det(J)# Jdet = J[0][0]*J[1][1]-J[0][1]*J[1][0]Jinv = np.linalg.inv(J)return Jinv,Jdetdef Bmatrix(r,s,Jinv):# 求BdNdr = diffNdr(r, s)dNds = diffNds(r, s)B1 = np.matrix([[1,0,0,0],[0,0,0,1],[0,<