3. MBB梁对称边界条件拓扑优化

MBB梁的对称性从何而来

MBB梁（Messerschmitt-Bölkow-Blohm梁）是拓扑优化领域的经典基准问题：矩形域上边界中部受集中力，左右下角简支。由于几何、载荷和支承都关于竖直中轴对称，最优材料分布也必然关于该轴对称。利用这一性质，只需对左半域建模，在中轴处施加对称边界条件，即可把设计变量数量减半。

对称边界条件的物理含义是：中轴上的点不允许有水平方向位移（u_x = 0），但允许竖直方向位移（u_y ≠ 0），因为梁在中轴处可以自由下挠。这与固定约束不同——固定约束会同时锁死两个方向，人为增加刚度、扭曲最优构型。

半模型建模要点

在有限元离散中，半模型的设计域为 [0, L/2] × [0, H]，右边界即对称轴。施加对称边界条件时，只需约束右边界所有节点的水平自由度 u_x = 0，竖直自由度保持自由。载荷施加在右上角节点，大小为全模型集中力的一半（F/2），因为对称轴只承担总载荷的一半。

体积约束需要按半模型重新计算：若全模型体积分数上限为 f，则半模型的材料体积上限为 f × (L/2 × H)，保持与全模型相同的材料用量比例。

可运行的Python实现

以下代码使用 scipy.sparse 组装刚度矩阵，以SIMP方法进行柔度最小化。设计域离散为 60×40 单元（半模型），体积分数上限 0.4。

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

# ---------- 参数 ----------
L, H = 60.0, 40.0          # 半模型尺寸 (mm)
nelx, nely = 60, 40        # 单元数
volfrac = 0.4              # 体积分数
penal = 3.0                # SIMP惩罚因子
E0, Emin = 1.0, 1e-9       # 材料模量
nu = 0.3                   # 泊松比
F_total = 1.0              # 全模型集中力
F_half = F_total / 2.0     # 半模型载荷

# ---------- 单元刚度矩阵 (Q4) ----------
def element_stiffness():
    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],
        [ 1/8+nu/8,  1/2-nu/6, -1/8+3*nu/8, -1/4-nu/12,
         -1/8-nu/8,   nu/6,    1/8-3*nu/8, -1/4+nu/12],
        [-1/4-nu/12, -1/8+3*nu/8, 1/2-nu/6,  1/8+nu/8,
           nu/6,      1/8-3*nu/8, -1/4+nu/12, -1/8-nu/8],
        [-1/8+3*nu/8, -1/4-nu/12, 1/8+nu/8,  1/2-nu/6,
          1/8-3*nu/8, -1/4+nu/12, -1/8-nu/8,   nu/6],
        [-1/4+nu/12, -1/8-nu/8,   nu/6,      1/8-3*nu/8,
          1/2-nu/6,  1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8],
        [-1/8-nu/8,   nu/6,      1/8-3*nu/8, -1/4+nu/12,
          1/8+nu/8,  1/2-nu/6, -1/8+3*nu/8, -1/4-nu/12],
        [  nu/6,      1/8-3*nu/8, -1/4+nu/12, -1/8-nu/8,
         -1/4-nu/12, -1/8+3*nu/8, 1/2-nu/6,  1/8+nu/8],
        [ 1/8-3*nu/8, -1/4+nu/12, -1/8-nu/8,   nu/6,
         -1/8+3*nu/8, -1/4-nu/12, 1/8+nu/8,  1/2-nu/6]
    ]) * E0 / (1 - nu**2)
    return k

KE = element_stiffness()

# ---------- 自由度编号 ----------
ndof = 2 * (nelx + 1) * (nely + 1)
def node_id(ix, iy):
    return iy * (nelx + 1) + ix

# ---------- 组装全局刚度矩阵 ----------
def assemble_K(x):
    rows, cols, vals = [], [], []
    for ey in range(nely):
        for ex in range(nelx):
            n1 = node_id(ex, ey)
            n2 = node_id(ex+1, ey)
            n3 = node_id(ex+1, ey+1)
            n4 = node_id(ex, ey+1)
            edof = [2*n1, 2*n1+1, 2*n2, 2*n2+1,
                    2*n3, 2*n3+1, 2*n4, 2*n4+1]
            ke = KE * (Emin + x[ey, ex]**penal * (E0 - Emin))
            for i in range(8):
                for j in range(8):
                    rows.append(edof[i])
                    cols.append(edof[j])
                    vals.append(ke[i, j])
    return coo_matrix((vals, (rows, cols)), shape=(ndof, ndof)).tocsr()

# ---------- 边界条件 ----------
def get_fixed_dofs():
    fixed = []
    # 左下角简支：约束 ux, uy
    n = node_id(0, 0)
    fixed += [2*n, 2*n+1]
    # 右下角简支（半模型右边界底部）：约束 uy
    n = node_id(nelx, 0)
    fixed += [2*n+1]
    # 对称轴（右边界）：约束 ux
    for iy in range(nely + 1):
        n = node_id(nelx, iy)
        fixed.append(2*n)
    return sorted(set(fixed))

fixed_dofs = get_fixed_dofs()
free_dofs = np.setdiff1d(np.arange(ndof), fixed_dofs)

# ---------- 载荷向量 ----------
F = np.zeros(ndof)
# 右上角节点施加 F/2 向下
n_load = node_id(nelx, nely)
F[2*n_load + 1] = -F_half

# ---------- SIMP主循环 ----------
x = np.full((nely, nelx), volfrac)
xPhys = x.copy()
loop, change = 0, 1.0
while change > 0.01 and loop < 200:
    loop += 1
    K = assemble_K(xPhys)
    U = np.zeros(ndof)
    U[free_dofs] = spsolve(K[free_dofs][:, free_dofs], F[free_dofs])
    # 柔度
    c = F @ U
    # 灵敏度
    dc = np.zeros((nely, nelx))
    for ey in range(nely):
        for ex in range(nelx):
            n1 = node_id(ex, ey)
            n2 = node_id(ex+1, ey)
            n3 = node_id(ex+1, ey+1)
            n4 = node_id(ex, ey+1)
            edof = [2*n1, 2*n1+1, 2*n2, 2*n2+1,
                    2*n3, 2*n3+1, 2*n4, 2*n4+1]
            Ue = U[edof]
            dc[ey, ex] = -penal * xPhys[ey, ex]**(penal-1) * (E0 - Emin) * (Ue @ KE @ Ue)
    # 灵敏度过滤（简化：3x3均值）
    dc_f = np.zeros_like(dc)
    for ey in range(nely):
        for ex in range(nelx):
            y0, y1 = max(0, ey-1), min(nely, ey+2)
            x0, x1 = max(0, ex-1), min(nelx, ex+2)
            dc_f[ey, ex] = dc[y0:y1, x0:x1].mean()
    # 最优准则更新
    l1, l2 = 1e-9, 1e9
    while (l2 - l1) / (l1 + l2) > 1e-3:
        lmid = 0.5 * (l1 + l2)
        xnew = np.maximum(0.001, np.maximum(x - 0.2,
                np.minimum(1.0, np.minimum(x + 0.2,
                x * np.sqrt(-dc_f / lmid)))))
        if xnew.mean() > volfrac:
            l1 = lmid
        else:
            l2 = lmid
    change = np.max(np.abs(xnew - x))
    x = xnew
    xPhys = x.copy()
    print(f"Iter {loop:3d}  c = {c:.4f}  vol = {xPhys.mean():.4f}  change = {change:.4f}")

print("优化完成，最终柔度 =", c)

结果验证与对称性检查

运行上述代码后，可观察以下关键指标：

• 柔度收敛曲线：柔度值应随迭代单调下降并趋于稳定，典型收敛在80~150次迭代之间。
• 体积分数：每步迭代后 xPhys.mean() 应稳定在 0.4 附近，偏差不超过 0.005。
• 构型对称性：将半模型结果沿对称轴镜像后，左右两侧材料分布应完全一致（数值误差范围内）。
• 计算效率：半模型自由度约为全模型的 1/2，刚度矩阵规模减小为 1/4，单次求解时间显著降低。