矩陣位移法——連續梁的結構剛度矩陣

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))