矩阵位移法——连续梁的结构刚度矩阵

import numpy as np
from math import sqrt, sin, cos, acos

np.set_printoptions(precision=3)
np.set_printoptions(suppress=True)  # 将科学数组转化为浮点数
EI = 1  # 为了方便计算，设EI=1

# n = input("请输入位移数量：")
n = 3
ElemNode = np.array([[0, 1], [1, 2], [2, 3]])  # 组成单元的两结点编号
NodeForce = np.array([[-60], [50], [30]])  # 作用在结点上的力，单位N·m
Ei = [0.75, 1.5, 1]
ki = np.array([[4 * EI, 2 * EI],
               [2 * EI, 4 * EI]])  # 单元刚度矩阵*L
NodeCoord = np.array([[0, 0], [6, 0], [14, 0], [20, 0]])  # 结点整体坐标

NodeX = NodeCoord[:, 0]  # 结点x坐标
NodeY = NodeCoord[:, 1]  # 结点y坐标


def k_i(i):
    """
    整体坐标系下单元刚度矩阵,i为单元编号
    :param i:结点编号
    :return: 整体坐标系下单元刚度矩阵，由于连续梁不用坐标变换，所以k = ke
    """
    NodeIndex = ElemNode[i, :]
    #print(NodeIndex)
    Node1_locx = NodeX[NodeIndex[0]]
    Node1_locy = NodeY[NodeIndex[0]]
    Node2_locx = NodeX[NodeIndex[1]]
    Node2_locy = NodeY[NodeIndex[1]]
    Len = sqrt((Node1_locx - Node2_locx) ** 2 + (Node1_locy - Node2_locy) ** 2)
   # print("len = %d " % (Len))
    return ki * Ei[i] / Len  # 假设为水平

def TransMatrix(angle):#坐标转换矩阵T
    return np.array([
        [ cos(angle),  sin(angle) , 0,      0       ,    0      ],
        [-sin(angle),  cos(angle),  0,      0       ,    0      ],
        [   0       ,      0     ,  1,      0       ,    0      ],
        [   0       ,      0     ,  0,    cos(angle), sin(angle)],
        [   0       ,      0     ,  0,   -sin(angle), cos(angle)],
        [   0       ,      0     ,  0,      0       ,    1      ],
])
# 对上述单元刚度矩阵的各元素，按照其行码和列码直接送入结构刚度矩阵，进行”对号入座“
NumberOfNodeFreeDof = 3
K = np.zeros((NumberOfNodeFreeDof, NumberOfNodeFreeDof))
for i in range(0, n):
    kii = k_i(i)
    # print(kii)
    for j in range(0, 2):
        for k in range(0, 2):
            if i + j - 1 < 0 or i + k - 1 < 0:
                continue
            K[i + j - 1, i + k - 1] += kii[j, k]
# print(K)

print("结构刚度矩阵K = ")
print(K)
print('-' * 20)
K_1 = np.linalg.inv(K)  # 求K矩阵的逆
print('K的逆 = ')
print(K_1)
print('-' * 20)
NodeDis = np.matmul(K_1, NodeForce)  # 🔺 = K-1*F
NodeDisplacement = np.zeros((1,n+1))

for i in range(0, n):
    NodeIndex = NodeDis[i,:]
    NodeDisplacement[0,i+1] = NodeIndex[0]
print(NodeDisplacement)
for i in range(0,n):
    kii = k_i(i)
    M1 = kii[0,0]*NodeDisplacement[0,i]+kii[0,1]*NodeDisplacement[0,i+1]
    M2 = kii[1, 0] * NodeDisplacement[0, i] + kii[1, 1] * NodeDisplacement[0, i + 1]
    print("在%d号杆件上，左端弯矩为M = %.3f, 右端弯矩为M = %.3f" % (i,M1,M2))