2D Quad4矩形单元计算实例编程[python]

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from collections import OrderedDict

# 一、定义函数
# 函数1：计算形函数矩阵N
def Q4_N(kesi, yita):
    N1 = (1-kesi)*(1-yita)/4
    N2 = (1+kesi)*(1-yita)/4
    N3 = (1+kesi)*(1+yita)/4
    N4 = (1-kesi)*(1+yita)/4
    N = np.array([
        [N1, 0, N2, 0, N3, 0, N4, 0],
        [0, N1, 0, N2, 0, N3, 0, N4]
    ])
    return N

# 函数2： 计算形矩阵G，用于计算应变矩阵
def Q4_G(kesi, yita):
    G = np.array([
        [yita/4 - 1/4, 0, 1/4 - yita/4, 0, yita/4 + 1/4, 0, - yita/4 - 1/4, 0],
        [kesi /4 - 1/4, 0, - kesi/4 - 1/4, 0, kesi/4 + 1/4, 0, 1/4 - kesi/4, 0],
        [0,yita/4 - 1/4, 0, 1/4 - yita/4, 0, yita/4 + 1/4, 0, - yita/4 - 1/4],
        [0, kesi /4 - 1/4, 0, - kesi/4 - 1/4, 0, kesi/4 + 1/4, 0, 1/4 - kesi/4]
    ])
    return G

# 函数3： 计算雅可比矩阵J
def Q4_J(kesi, yita, xy):
    NN = np.array([
        [-(1-yita)/4, (1-yita)/4, (1+yita)/4, -(1+yita)/4],
        [-(1-kesi)/4, -(1+kesi)/4, (1+kesi)/4, (1-kesi)/4]
    ])
    J = NN @ xy
    return J

# 二、输入信息
# 定义结构材料、截面属性
E = 2.000E+08   # 弹性模量
v = 0.3     # 泊松比
L = 2       # 长度
H = 0.5     # 高度
t = 0.2    # 宽度

# 网格划分
DoE = 0.25  # 网格尺寸
NEofH = int(H/DoE)  # 高度方向单元数
NEofL = int(L/DoE)  # 长度方向单元数
NodeCoord = np.zeros(((NEofH+1)*(NEofL+1), 2))  # 节点坐标矩阵初始化
EleNode = np.zeros([int(NEofH*NEofL), 4])      # 单元节点矩阵初始化

# 定义节点坐标
for i in range(1, int(NEofL)+2):
    for j in range(1, int(NEofH)+2):
        NodeCoord[((i-1)*(int(NEofH)+1)+j)-1, 0] = (L/NEofL)*(i-1)
        NodeCoord[((i-1)*(int(NEofH)+1)+j)-1, 1] = (H/NEofH)*(j-1)

# 定义单元节点
for i in range(1, int(NEofL)+1):
    for j in range(1, int(NEofH)+1):
        EleNode[((i-1)*int(NEofH)+j)-1, 0] = (NEofH+1)*(i-1) + 1+(j-1)
        EleNode[((i-1)*int(NEofH)+j)-1, 1] = (NEofH+1)*i + 1+(j-1)
        EleNode[((i-1)*int(NEofH)+j)-1, 2] = (NEofH+1)*i + 2+(j-1)
        EleNode[((i-1)*int(NEofH)+j)-1, 3] = (NEofH+1)*(i-1) + 2+(j-1)

# 单元和节点数量
numEle = EleNode.shape[0]  # 单元数量
numNode = NodeCoord.shape[0]  # 节点数量

# 节点坐标值向量
xx = NodeCoord[:, 0]
yy = NodeCoord[:, 1]

# 自由度数量
numDOF = 2 * numNode    
# 约束自由度
ConstrainedDOF = np.zeros([1, 2*(int(NEofH)+1)], dtype='int')
for i in range(1, 2*(int(NEofH)+1)+1):
    ConstrainedDOF[0][i -1] =i

#  定义整体位移向量
displacement = np.zeros((1, numDOF))

# 定义整体载荷向量
force = np.zeros((1, numDOF))
p = -100
force[0][int(numDOF)-int(NEofH)-1] = p

# 积分点及积分点权重
kesi_yita = np.array([
    [-1/np.sqrt(3), -1/np.sqrt(3)],
    [1/np.sqrt(3), -1/np.sqrt(3)],
    [1/np.sqrt(3), 1/np.sqrt(3)],
    [-1/np.sqrt(3), 1/np.sqrt(3)]
])
w = np.array([1, 1])
ww = np.array([w[0]*w[0], w[1]*w[0], w[1]*w[1], w[0]*w[1]])

# 计算分析
# 定义总体刚度矩阵
stiffness_sum = np.zeros((numDOF, numDOF))

# 定义弹性矩阵D
D = (E/(1-v**2))*np.array([
    [1, v, 0],
    [v, 1, 0],
    [0, 0, (1-v)/2]
])


# 遍历每个单元，计算单元刚度矩阵并组装到总体刚度矩阵中
for i in range(numEle):
    noindex = EleNode[i]
    xy = np.zeros([4,2])
    for j in range(4):
        xy[j, 0] = xx[int(noindex[j])-1]
        xy[j, 1] = yy[int(noindex[j])-1]
    eleK = np.zeros((8, 8))
    for m in range(4):
        kesi = kesi_yita[m][0]
        yita = kesi_yita[m][1]
        J = Q4_J(kesi, yita, xy)
        A = (1/np.linalg.det(J))*np.array([
            [J[1,1], -1*J[0,1], 0, 0],
            [0, 0, -1*J[1,0], J[0,0]],
            [-1*J[0,1], J[0,0], J[1,1], -1*J[1,0]]
        ])
        G = Q4_G(kesi, yita)
        B = A @ G
        eleK =eleK+ww[m]*B.T @ D @ B * t * np.linalg.det(J)

    # 根据单元节点编号，确定单元自由度
    eleDof = np.array([
        noindex[0]*2-2, noindex[0]*2-1,
        noindex[1]*2-2, noindex[1]*2-1,
        noindex[2]*2-2, noindex[2]*2-1,
        noindex[3]*2-2, noindex[3]*2-1
    ])
    # 将单元刚度矩阵组装到总体刚度矩阵中
    for i in range(8):
        for j in range(8):
            stiffness_sum[int(eleDof[i]), int(eleDof[j])] += eleK[i][j]

# 根据给定的约束条件，删除指定被约束自由度对应的行和列
stiffness_sum_sim = np.delete(stiffness_sum, ConstrainedDOF-1, 0)
stiffness_sum_sim = np.delete(stiffness_sum_sim, ConstrainedDOF-1, 1)

# 解方程，求位移
activeDof = np.delete(np.arange(numDOF) +1, ConstrainedDOF-1)
displacement[0][activeDof -1] = np.array(np.asmatrix(stiffness_sum_sim).I * np.asmatrix(force[0][activeDof -1]).T).T

# 按节点编号输出两个方向的节点位移分量，格式：节点号，位移1，位移2
Nodeindex = np.arange(1, numNode+1).reshape(numNode, 1)
tempdisp =np.hstack((Nodeindex, displacement.reshape(numNode, 2)))

# 计算节点应力
stress_Node = np.zeros((numEle * 4, 5))
kesi_yita_Node = np.array([
    [-1, -1],
    [1, -1],
    [1, 1],
    [-1, 1] 
])
for i in range(numEle):
    noindex = EleNode[i]
    xy_node = np.zeros((4,2))
    d = np.zeros((8,1))

    for j in range(4):
        xy_node[(j,0)] = xx[int(noindex[j])-1]
        xy_node[(j,1)] = yy[int(noindex[j])-1]
        d[int(j*2)] = displacement[0][int(noindex[j]-1)*2]
        d[int(j*2+1)] = displacement[0][int(noindex[j]-1)*2+1]
    for m in range(4):
        kesi = kesi_yita_Node[m][0]
        yita = kesi_yita_Node[m][1]
        J = Q4_J(kesi, yita, xy)
        A = (1/np.linalg.det(J))*np.array([
            [J[1,1], -1*J[0,1], 0, 0],
            [0, 0, -1*J[1,0], J[0,0]],
            [-1*J[1,0], J[0,0], J[1,1], -1*J[0,1]]
        ])
        G = Q4_G(kesi, yita)
        B = A @ G
        sima = D @ B @ d
        stress_Node[i*4+m] = np.asarray([i+1, noindex[m], sima[0][0], sima[1][0], sima[2][0]], dtype=object)

# 计算支座反力
reation =np.zeros((ConstrainedDOF.size, 1))
for i in range(ConstrainedDOF.size):
    rownum = ConstrainedDOF[0][i]
    reation[i] = stiffness_sum[rownum -1,] @ displacement.T - force[0][rownum -1]


# 四、计算用于绘图相关参数
# 计算个单元最大长度maxlen
maxlen = -1
for i in range(numEle):
    noindex = EleNode[i]
    deltax = xx[int(noindex[1]-1)] - xx[int(noindex[0]-1)]
    deltay = yy[int(noindex[1]-1)] - yy[int(noindex[0]-1)]
    L = np.sqrt(deltax**2 + deltay**2)
    if L > maxlen:
        maxlen = L

# 计算节点最大绝对位移
maxabsDisp = np.sqrt((displacement[0][0] * displacement[0][0])+(displacement[0][1] * displacement[0][1]))
for i in range(numNode):
    indexU = 2 * (i+1) -2
    indexV = 2 * (i+1) -1
    dispU = displacement[0][indexU]
    dispV = displacement[0][indexV]
    absDisp = np.sqrt(dispU**2 + dispV**2)
    if absDisp > maxabsDisp:
        maxabsDisp = absDisp

# 计算变形图变形放大系数
factor = 0.1
scalefactor = 0.2
if maxabsDisp > 1e-30:
    factor = scalefactor * maxlen / maxabsDisp

# 五、绘图
fig = plt.figure(figsize=(15, 6))
plt.subplots_adjust(left=0.1, right = 1)

# 绘制为变形的结构线条
for i in range(numEle):
    noindex = np.concatenate([ [EleNode[i][0]], [EleNode[i][1]], [EleNode[i][2]], [EleNode[i][3]],  [EleNode[i][0]],])
    x = xx[noindex[:].astype('int')-1]
    y = yy[noindex[:].astype('int')-1]
    line1 = plt.plot(x, y, label='Undeformed Shape', color='grey', linestyle='--')

# 绘制变形后的结构
xx_Deformed = np.empty_like(xx)
yy_Deformed = np.empty_like(yy)

for i in range(numEle):
    #获取第i个单元的节点索引，闭合曲线需要重复第一个节点方向索引
    noindex = np.concatenate([[EleNode[i][0]], [EleNode[i][1]], [EleNode[i][2]], [EleNode[i][3]], [EleNode[i][0]]])
    indexU= 2*noindex[:].astype('int')-2
    indexV= 2*noindex[:].astype('int')-1 
    dispU= displacement[0][indexU]
    dispV= displacement[0][indexV]
    xxDeformed = xx[noindex[:].astype('int')-1] + dispU*factor
    yyDeformed = yy[noindex[:].astype('int')-1] + dispV*factor
    xx_Deformed[noindex.astype('int')-1] = xxDeformed
    yy_Deformed[noindex.astype('int')-1] = yyDeformed
    middlex = 0.5*(xxDeformed[1]+xxDeformed[0])
    middley = 0.5*(yyDeformed[1]+yyDeformed[0])
    line2 = plt.plot(xxDeformed, yyDeformed, label='Deformed Shape', color='black', marker='o', mfc='w')


#标记节点编号
for i in range(numNode):
    plt.text(xx_Deformed[i]+0.03, yy_Deformed[i], str(i+1), color='red', ha='center', va='center')


#设置坐标格式，使得图形以适合的比例进行显示
max_xx=max(max(xx), max(xx_Deformed))
max_yy=max(max(yy), max(yy_Deformed))
min_xx=min(min(xx),min(xx_Deformed))
min_yy=min(min(yy),min(yy_Deformed))

minx= -0.2*(max_xx-min_xx)+min_xx
maxx=0.2*(max_xx-min_xx)+max_xx
miny= -0.2*(max_yy-min_yy)+min_yy
maxy=0.2*(max_yy-min_yy)+max_yy

plt.xlim(0.2*minx,0.89*maxx)
plt.ylim(miny, maxy)
#设置x坐标轴刻度
plt.xticks(np.arange(0, 2.2,0.2))
#设置x坐标轴标签
plt.xlabel('X')
#设置y坐标轴标签
plt.ylabel('Y')


#绘制应力云图
sima_x= np.zeros((numNode, 1))
for stress_Node_sl in stress_Node:
    sima_x[int(stress_Node_sl[1])-1]= stress_Node_sl[2]#将方向应力数组赋值给变量2
Z=sima_x
#将x坐标数组展平为一维数组#将y坐标数组展平为一维数组
x=xx_Deformed.flatten()
y=yy_Deformed.flatten()
#将应力数组展平为一维数组
z=Z.flatten()
#构建三角剖分对象
triang =tri.Triangulation(x, y)
tt=plt.tricontourf(triang,z,cmap="viridis")
plt.colorbar(tt)


#去掉重复标签
#绘制应力云图的颜色条
handles, labels = plt.gca().get_legend_handles_labels()#获取图例句柄和标签
by_label =OrderedDict(zip(labels,handles))#按照标签和句柄创建有序字典
plt.legend(by_label.values(),by_label.keys(), loc='upper right')#绘制图例，并设置位置为右上角
plt.show()