2. 简支梁体积约束下的刚度最大化

简支梁刚度最大化：从工况到最优构型

本例对象为跨度 L=120 mm、高 H=40 mm 的二维简支梁，左端下缘铰支（约束 x、y 位移），右端下缘滚动支承（约束 y 位移），跨中上缘作用集中力 F=100 N，方向竖直向下。设计域离散为 120×40 的四边形单元网格，材料弹性模量 E=1（归一化），泊松比 ν=0.3，体积分数上限 f 分别取 0.3、0.4、0.5 进行对比。

刚度最大化等价于柔度最小化。采用 SIMP 插值模型，单元刚度按惩罚因子 p=3 在实体与空洞之间过渡，优化列式为：

min  c(ρ) = U^T K U = Σ_e (ρ_e)^p · u_e^T k0 u_e
s.t. K U = F
     Σ_e ρ_e · v_e / V0 ≤ f
     0 < ρ_min ≤ ρ_e ≤ 1

其中 ρ_e 为单元相对密度，ρ_min=0.001 防止刚度矩阵奇异，v_e 为单元体积，V0 为设计域总体积。灵敏度由伴随法得到：∂c/∂ρ_e = −p·(ρ_e)^(p−1)·u_e^T k0 u_e。采用最优准则（OC）更新密度，并配合灵敏度过滤（过滤半径 r_min=2 个单元）抑制棋盘格。

实现要点与结果判读

在 MATLAB 或 Python 中实现时，核心循环为：组装全局刚度矩阵 → 施加边界条件 → 求解位移 → 计算柔度与灵敏度 → 灵敏度过滤 → OC 更新密度 → 检查收敛。以下为关键代码骨架（Python 风格，使用 NumPy 与 SciPy 稀疏求解）：

import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve

# 网格与材料参数
nelx, nely = 120, 40
E0, nu, p, rmin = 1.0, 0.3, 3.0, 2.0
volfrac = 0.4
penal = 3.0

# 单元刚度矩阵（平面应力，单位厚度）
def lk():
    k = np.array([1/2-nu/6, 1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8,
                  -1/4+nu/12, -1/8-nu/8, nu/6, 1/8-3*nu/8])
    KE = 1/(1-nu**2)*np.array([
        [k[0],k[1],k[2],k[3],k[4],k[5],k[6],k[7]],
        [k[1],k[0],k[7],k[6],k[5],k[4],k[3],k[2]],
        [k[2],k[7],k[0],k[1],k[6],k[3],k[4],k[5]],
        [k[3],k[6],k[1],k[0],k[7],k[2],k[5],k[4]],
        [k[4],k[5],k[6],k[7],k[0],k[1],k[2],k[3]],
        [k[5],k[4],k[3],k[2],k[1],k[0],k[7],k[6]],
        [k[6],k[3],k[4],k[5],k[2],k[7],k[0],k[1]],
        [k[7],k[2],k[5],k[4],k[3],k[6],k[1],k[0]]])
    return KE

# 初始化密度与过滤
x = np.ones(nely*nelx)*volfrac
xmin, xmax = 0.001, 1.0
loop, change = 0, 1.0
while change > 0.01 and loop < 200:
    loop += 1
    # 组装刚度矩阵（SIMP插值）
    # ... 此处省略组装细节，实际需按单元循环组装 K
    # 求解位移 U
    # 计算柔度 c 与灵敏度 dc
    # 灵敏度过滤
    # OC 更新 xnew
    change = np.max(np.abs(xnew - x))
    x = xnew
    print(f"Iter {loop}: c = {c:.4f}, vol = {np.mean(x):.4f}")

不同体积分数的对比结论

从迭代曲线看，柔度在前期快速下降，约 40~60 步后趋于平稳；体积分数越大，收敛后的柔度越低，但材料利用率（单位体积刚度增益）并非线性提升。工程中需在刚度需求与减重目标之间折中选取 f。