二维4单元9节点总刚度组装以及9个节点的力和位移求解-节点编号对结果的影响

本人参照力学猿https://blog.csdn.net/Youngist/article/details/106675098的python编程，在只改变节点编号的情况下（节点编号更改如下图），对源程序进行以下修改

1. 将func_k函数对应的x0,y0进行更改

def func_k(i, j, k, l):  # 根据参数得到4个单元刚度矩阵
    intPoint = [-1 / sqrt(3), 1 / sqrt(3)]  # 二阶高斯点坐标
    weight = [1.0, 1.0]  # 权重
    E = 1  # 弹性模量
    nu = 0.3  # 泊松比
    C = planeStressC(E, nu)
    x0 = [-1, 1, 1, -1, -1, 0, 1, 0, 0]
    y0 = [-1, -1, 1, 1, 0, -1, 0, 1, 0]
    x = [x0[i - 1], x0[j - 1], x0[k - 1], x0[l - 1]]
    y = [y0[i - 1], y0[j - 1], y0[k - 1], y0[l - 1]]
    K = np.zeros((8, 8))
    for j in range(0, 2):  # 数值积分求K
        for i in range(0, 2):
            r = intPoint[i]
            s = intPoint[j]
            Jinv, Jdet = jacobian(x, y, r, s)
            B = Bmatrix(r, s, Jinv)
            BT = np.transpose(B)
            tmp = BT * C * B * Jdet
            K = K + tmp
    return K

2.更改K_element对象对应的不同i,jk,l对应的实参

SE_1 = K_element(1, 6, 9, 5)
print("单元1的刚度矩阵：{}".format(SE_1.K))
SE_2 = K_element(6, 2, 7, 9)
print("单元2的刚度矩阵：{}".format(SE_2.K))
SE_3 = K_element(9, 7, 3, 8)
print("单元3的刚度矩阵：{}".format(SE_3.K))
SE_4 = K_element(5, 9, 8, 4)
print("单元4的刚度矩阵：{}".format(SE_4.K))

3.改变对应求节点力和结点位移的对应程序

new2_F = np.matrix([[1000], [0], [-1000], [0], [-1000], [0], [0], [-1000], [1000], [0], [0], [-1000], [0], [-1000]])
new_kk_row = np.delete(kk, [0, 1, 2, 3], axis=0)
new_kk = np.delete(new_kk_row, [0, 1, 2, 3], axis=1)  # 在原总刚度下删除U=0所在的行和列，得到一个新的矩阵
new2_U = my_LUsolve(new_kk, new2_F)
U1 = np.zeros((18, 1))
for i in range(14):
    U1[i + 4, 0] = new2_U[i, 0]
print("各节点的位移U1为：")
print(U1)
print(np.shape(kk))
print(np.shape(U1))
print("各节点的力F1为：")
F1 = kk * np.matrix(U1)
print(F1)

其它的程序和源程序一样，输出结果如下

源程序输出结果如下

可以明显看到输出的节点力和结点位移不一样，在固定支撑的节点力不一样可以用应力奇异点来解释，但对于实际位置的结点位移应该不会随节点编号的改变而改变改变，比如原本编号对应的1，2、3节点位移，应该分别和改变编号后的1，6，2节点位移一样，那么造成结果不一样的原因是什么，能帮本人解答一下吗