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
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,N4
def diffNdr(r,s):
dN1dr = 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 dNdr
def diffNds(r,s):
dN1ds = 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 dNds
def jacobian(x,y,r,s):
dNdr = 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)
Jinv = np.linalg.inv(J)
return Jinv,Jdet
def Bmatrix(r,s,Jinv):
dNdr = diffNdr(r, s)
dNds = diffNds(r, s)
B1 = np.matrix([[1,0,0,0],[0,0,0,1],[0,<