Python有限元代码组装整体刚度矩阵时出现广播错误求助

错误原因 #

报错核心是节点索引逻辑错误：代码误用1-based节点编号规则，但Python数组采用0-based索引。当i=0（第一个单元）时，node1=(0-1)*2+1=-1，导致K_global[node1:node3+1, node1:node3+1]切片结果为空矩阵（形状(0,0)），无法与3×3的单元刚度矩阵执行加法运算，触发广播不兼容的ValueError。

同时代码还存在其他索引问题：

• 单元刚度矩阵循环中，错误用i-1索引单元面积和长度，导致第一个单元取到最后一个元素的值
• 载荷组装循环的节点索引沿用错误的1-based逻辑
• 应力计算循环跳过第一个单元，导致结果不完整

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

E = 2e10  # Young's Modulus (kg/m^2)
L = 28    # Length of the bar (m)
while True:
    # Define Problem
    no_elements = int(input("Enter the number of elements: "))
    no_nodes_per_element = 3
    no_nodes = no_elements * (no_nodes_per_element - 1) + 1
    print(no_nodes)
    element_length = L / no_elements
    element_length = np.ones(no_elements) * element_length
    mid = np.cumsum(element_length) - element_length / 2
    node_length = L / (no_nodes-1)
    A = lambda x: (11/19600) * x**2 - (19/700) * x + 0.36  # Area equation
    elem_areas = A(mid)

    # Stiffness Matrix of elements
    elementStiffnessMatrices = []
    for i in range(no_elements):
        A_e = elem_areas[i]  # 修正索引为i，而非i-1
        E_e = E
        Le = element_length[i]  # 修正索引为i，而非i-1
        # Modified stiffness matrix for a truss element with 3 nodes
        Ke_element = (A_e * E_e / (3 * Le)) * np.array([[7, -8, 1], [-8, 16, -8], [1, -8, 7]])
        elementStiffnessMatrices.append(Ke_element)

    # Global Stiffness Matrix
    K_global = np.zeros((no_nodes, no_nodes))
    for i in range(no_elements):
        # 修正为0-based节点索引
        node1 = i * 2
        node2 = i * 2 + 1
        node3 = i * 2 + 2
        K_global[node1:node3+1, node1:node3+1] += elementStiffnessMatrices[i]

    K_original = K_global.copy()
    print(K_original)

    # Forces
    appliedLoadpermeter = 2000  # Given line load is 2000 kg/m
    EEnodalForces = np.zeros((no_elements, 3))

    for i in range(no_elements):
        Le = element_length[i]
        # Energy Equivalent Nodal Forces for truss element with 3 nodes
        EEnodalForces[i, :] = (appliedLoadpermeter * Le / 2) * np.array([1, 1, 1])
    
    globalEENodalForces = np.zeros(no_nodes)
    for i in range(no_elements):
        # 修正为0-based节点索引
        node1 = i * 2
        node2 = i * 2 + 1
        node3 = i * 2 + 2
        indices = [node1, node2, node3]
        globalEENodalForces[indices] += EEnodalForces[i, :]
    
    F_original = globalEENodalForces.copy()

    # Applying Boundary conditions
    Fixedsupport = 1
    # Delete the rows and columns corresponding to the fixed support
    K_global = np.delete(K_global, 0, axis=0)
    K_global = np.delete(K_global, 0, axis=1)
    globalEENodalForces = np.delete(globalEENodalForces, 0)

    # Solve for displacements
    Displacement = np.linalg.solve(K_global, globalEENodalForces)
    Displacement = np.insert(Displacement, 0, 0)

    # Stress
    elem_Stress = np.zeros((no_elements, 3))
    for i in range(no_elements):  # 修正循环范围，包含所有单元
        u1 = Displacement[i * 2]
        u2 = Displacement[i * 2 + 1]
        u3 = Displacement[i * 2 + 2]
        Le = element_length[i]
        Stress = E * np.array([[-3/Le, 4/Le, -1/Le],
                               [-1/Le, 0, 1/Le],
                               [1/Le, -4/Le, 3/Le]]) @ np.array([u1, u2, u3])
        elem_Stress[i, 0] = Stress[0]
        elem_Stress[i, 1] = Stress[1]
        elem_Stress[i, 2] = Stress[2]

################################################################################
# Step-9: Exact Values
# Equation of the displacment for Exact solution
    def U_x(x):
        return - (49 * np.log(11 * x**2 - 532 * x + 7056)) / 550000 - (21 * 5**(1/2) * 7**(1/2) * np.arctan((19 * 5**(1/2) * 7**(1/2)) / 35 - (11 * 5**(1/2) * 7**(1/2) * x) / 490)) / 1375000 + 9.0415e-4

#Equation of the Stress for Exact solution
    def S_x(x) : return (-2000*x + 56000)/((11/19600)*x**2 - (19/700)*x + 0.36)

################################################################################
# Displacement figure
    AXIS = np.arange(0, L + node_length, node_length)
    plt.figure()
    plt.plot(AXIS, Displacement, "b")
    plt.plot(np.arange(0, L + 0.1, 0.1), U_x(np.arange(0, L + 0.1, 0.1)), "r")
    plt.legend(["FE", "Exact"])
    plt.xlabel("Length of element (m)")
    plt.ylabel("Displacement (m)")
    plt.xlim([0, 28])
    plt.show()

# Stress figure
    plt.figure(2)
    for i in range(no_elements):
        a = mid[i] - element_length[i] / 2
        b = mid[i] + element_length[i] / 2
        plt.plot(np.linspace(a, b, 100), elem_Stress[i, 1] * np.ones(100), "b", label="FE" if i==0 else "", linewidth=2)   
    plt.plot(np.arange(0, L + 0.1, 0.1), S_x(np.arange(0, L + 0.1, 0.1)), "r", label="Exact", linewidth=2)
    plt.legend()
    plt.xlabel("Length of the element (m)")
    plt.ylabel("Stress (kg/m2)")
    plt.title("Comparison of FE Stress and Exact Solution")
    plt.show()

关键修正点说明 #

1. 单元刚度矩阵索引修正：将elem_areas[i-1]和element_length[i-1]改为elem_areas[i]和element_length[i]，匹配0-based数组索引
2. 全局刚度矩阵组装索引修正：将节点计算逻辑改为node1=i*2、node2=i*2+1、node3=i*2+2，确保每个单元的3个节点对应全局矩阵的正确切片范围
3. 载荷组装索引修正：同步修正载荷组装循环中的节点索引，确保载荷正确分配到对应节点
4. 应力计算循环修正：将循环范围从range(1, no_elements)改为range(no_elements)，并修正位移索引为0-based，确保所有单元的应力都被计算
5. 图例重复标注修正：在应力绘图循环中，仅在第一个单元添加"FE"图例，避免重复生成图例

内容的提问来源于stack exchange，提问作者user22989254