3. 平面梁单元弯曲问题有限元实现

欧拉-伯努利梁单元的核心假设与刚度矩阵

平面梁单元与杆单元的本质区别在于：杆单元只有轴向自由度，而梁单元每个节点包含横向位移 w和截面转角 θ两个自由度。欧拉-伯努利梁理论建立在以下假设之上：

• 变形前垂直于中面的横截面，变形后仍保持平面且垂直于变形后的中面（直法线假设）；
• 忽略剪切变形的影响，即截面转角等于挠度的一阶导数：θ = dw/dx；
• 梁为等截面直梁，材料在线弹性范围内工作。

基于上述假设，两节点梁单元（长度 L，抗弯刚度 EI）的节点位移向量为 d = [w₁, θ₁, w₂, θ₂]ᵀ，对应的单元刚度矩阵为：

k = (EI/L³) ×
[  12      6L     -12      6L  ]
[  6L     4L²     -6L     2L²  ]
[ -12     -6L      12     -6L  ]
[  6L     2L²     -6L     4L²  ]

该矩阵的物理含义：对角线元素代表单位位移引起的节点力，非对角线元素代表耦合效应。例如 k₁₁ = 12EI/L³ 表示节点1产生单位横向位移（θ₁=0, w₂=0, θ₂=0）时所需的横向力；k₁₂ = 6EI/L² 表示节点1单位转角引起的节点1横向力。

有限元求解流程与关键实现

以悬臂梁为例，取 4 个单元（5 个节点），每个节点 2 个自由度，总自由度数为 10。整体刚度矩阵由各单元刚度矩阵按节点自由度编号组装而成。施加边界条件时，固定端节点 1 的 w₁ 和 θ₁ 被约束，对应行列从方程组中消去。

求解线性方程组 K·d = F 后，可提取自由端节点的挠度和转角。与解析解对比时需注意：有限元节点位移是精确的（对于欧拉-伯努利梁单元，节点位移在无分布载荷时与解析解完全一致），但单元内部的挠度需要形函数插值，存在一定误差。

Python 实现与结果验证

以下代码实现悬臂梁的有限元求解，并与解析解对比。梁参数：L=1m，EI=1000 N·m²，端部集中力 P=10 N。

import numpy as np

def beam_element_stiffness(EI, L):
    """返回两节点梁单元的4x4刚度矩阵"""
    k = EI / L**3 * np.array([
        [ 12,   6*L,  -12,   6*L],
        [ 6*L,  4*L**2, -6*L, 2*L**2],
        [-12,  -6*L,   12,  -6*L],
        [ 6*L,  2*L**2, -6*L, 4*L**2]
    ])
    return k

def solve_cantilever(L_total, EI, P, n_elem):
    """悬臂梁有限元求解，返回节点位移向量"""
    L = L_total / n_elem
    n_node = n_elem + 1
    n_dof = 2 * n_node
    K = np.zeros((n_dof, n_dof))
    k_e = beam_element_stiffness(EI, L)

    # 组装整体刚度矩阵
    for e in range(n_elem):
        dofs = [2*e, 2*e+1, 2*e+2, 2*e+3]
        for i in range(4):
            for j in range(4):
                K[dofs[i], dofs[j]] += k_e[i, j]

    # 载荷向量：端部集中力 P 作用于最后一个节点的横向自由度
    F = np.zeros(n_dof)
    F[-2] = -P  # 向下为负

    # 边界条件：节点1的 w 和 theta 固定
    fixed = [0, 1]
    free = [i for i in range(n_dof) if i not in fixed]

    # 求解
    d = np.zeros(n_dof)
    d[free] = np.linalg.solve(K[np.ix_(free, free)], F[free])
    return d, n_node

# 参数
L_total = 1.0
EI = 1000.0
P = 10.0
n_elem = 4

d, n_node = solve_cantilever(L_total, EI, P, n_elem)

# 提取自由端挠度和转角
w_tip = d[-2]
theta_tip = d[-1]

# 解析解
w_exact = -P * L_total**3 / (3 * EI)
theta_exact = -P * L_total**2 / (2 * EI)

print(f"单元数: {n_elem}, 节点数: {n_node}")
print(f"自由端挠度  FEM: {w_tip:.6e} m,  解析: {w_exact:.6e} m")
print(f"自由端转角  FEM: {theta_tip:.6e} rad, 解析: {theta_exact:.6e} rad")
print(f"挠度相对误差: {abs((w_tip - w_exact)/w_exact)*100:.4f}%")
print(f"转角相对误差: {abs((theta_tip - theta_exact)/theta_exact)*100:.4f}%")

运行上述代码，预期输出如下（节点位移在集中力作用下与解析解一致，误差来自浮点运算）：

单元数: 4, 节点数: 5
自由端挠度  FEM: -3.333333e-03 m,  解析: -3.333333e-03 m
自由端转角  FEM: -5.000000e-03 rad, 解析: -5.000000e-03 rad
挠度相对误差: 0.0000%
转角相对误差: 0.0000%

简支梁均布载荷的等效节点载荷

对于长度为 L 的梁单元，承受均布载荷 q 时，等效节点载荷向量为：

f_eq = [qL/2, qL²/12, qL/2, -qL²/12]ᵀ

其中前两个分量对应节点 1 的横向力和弯矩，后两个对应节点 2。将各单元的等效节点载荷组装到整体载荷向量后，再施加简支边界条件（两端横向位移为零），即可求解。