Page 238 - 《振动工程学报》2025年第11期

2696 振 动 工 程 学 报 第 38 卷

投影 [35] ： 2.2 拓扑优化方法
∑
定义目标频段内的超材料的平均反射系数为
H i,j ξ i

j∈N d
ξ i = ，H i,j = max(0,r −
x i − x j

) （32） R = mean(R( f))，为使得优化结果的吸声系数尽可能

∑
H i,j
大，目标函数设置为最小化平均反射系数。由于鲁
j∈N d
棒性公式会使得优化过程中产生三个不同的设计实
tanh(βη)+tanh(β(ξ i −η))
ˆ ξ i = （33） 现。因此，该优化问题的目标函数可转换为关于设
tanh(βη)+tanh(β(1−η))
式中， x i 表示第 i 个单元网格的中心位置； H i,j 表示根 计 变量 ξ 和 旋 转 角 θ 的 最 大 值 最 小 化 问 题 （ min-
max），即
据指定的过滤半径 r 决定的权重函数；η 和 β 分别表
min max(R e ,R b ,R d ) （36）
示投影阈值和投影强度。 ξ,θ
通过式 (32) 的密度过滤后获得的中间变量 ξ i 通 将该最大值最小化问题转换为约束公式，同时
常为灰度单元，即其值在 0~1 之间。为推动设计变 在拓扑优化过程中加入刚度约束和密度约束，可以
量转变为 0 或 1 以形成明确的物理结构，需要通过 得到如下的优化公式：

 min
式 (33) 中 的 Heaviside 投 影 操 作 进 一 步 将 中 间 变 量  ξ,θ Ψ






ξ i 投影为物理变量 。  s.t. R e −Ψ ⩽ 0，R b −Ψ ⩽ 0，R d −Ψ ⩽ 0

ˆ ξ i


 /E t +1 ⩽ 0

−E x, ˆ ξ b （37）
为保证优化过程的稳健性，本文优化过程中使 


 /E t +1 ⩽ 0
 −E y, ˆ ξ b

用了鲁棒性公式，即  /V t −1 ⩽ 0



 Vˆ ξ b


 0 ⩽ Ψ ⩽ 1，0 ⩽ ξ ⩽ 1，0 ⩽ θ ⩽ π/2
式(31) 式(32) 式(33)
ˆ ˆ ˆ
ξ −−−−→ ξ −−−−→ ξ −−−−−−−→ (ξ e ,ξ b ,ξ d ) （34）
L r β,(η e ,η b ,η d ) 式 中， E x, ˆ ξ b 表 示 蓝 图 设 计 实 现 的 x 方 向 等 效 刚 度 ；
ˆ
ˆ
ˆ
其中， ξ e 、ξ b 、ξ d 分别表示侵蚀设计、蓝图设计和膨 E y, ˆ ξ b 表示 y 方向等效刚度； E t 为指定的刚度约束阈
，
胀 设 计； (η e ,η b ,η d ) = (η b +∆η,η b ,η b −∆η) ∆η为 侵 蚀 设 值，通过结构刚度约束，可以确保实现的结构在静水
计或者膨胀设计与蓝图设计的投影阈值差值。 压下具有小变形量，从而展现出良好的耐压性能；
侵蚀设计倾向于“消除”小尺寸特征，使结构更 Vˆ ξ b 表示蓝图设计实现的结构体积分数； V t 为指定的
保守、材料分布更集中，从而抑制数值不稳定；膨胀 体积分数上限，通过体积分数约束可确保较小的结
设计倾向于“扩展”材料区域，使结构更健壮、对不 构等效密度，实现轻质的超材料。由于反射系数取
确定性更不敏感；蓝图设计则介于两者之间，反映优 值范围为 [0 1]，目标值 Ψ 取值范围与其相同。通过
化过程中得到的理论最优拓扑结构形态，因此最终 该目标值可以有效约束三个不同设计的反射系数值。

元胞形态结构使用蓝图设计。该过程元胞结构投影
2.3 逆向设计过程
示意图如图 7 所示。通过鲁棒性公式，可以提高拓
扑优化设计过程中的结构稳定性。 为设计具有高刚度以及轻质特性的宽频带吸声
基于各向同性材料罚函数模型，每个网格单元 超材料，本文开展了结构拓扑优化设计，流程图如图 8
内的杨氏模量 E 和密度 ρ 描述为物理变量 ˆ ξ i 的插值 所示。
函数，即 首先，预设优化参数，元胞结构离散为 N d =100×
 p 1 100 个四边形网格单元。式 (34) 中，r=0.08a，η e =0.55，
 ρ( ˆ ξ i ) = ρ 0 + ˆ ξ (ρ 1 −ρ 0 )

 i
 （35） η b =0.50，η d =0.45；式 (35) 中，惩罚因子 p 1 =1，p 2 =3.5。约


 p 2
 E( ˆ ξ i ) = E 0 + ˆ ξ (E 1 − E 0 )
i 束阈值中，设置 E t =150 MPa，V t =0.55。随后，元胞结
式中，ρ 1 和 E 1 分别为基体材料的杨氏模量和密度； 构设置为正交各向异性，将设计域基于式 (34) 进行
ρ 0 = 10 ρ 1 和 E 0 = 10 E 1 为人工变量参数；p 1 和 p 2 分 处理为物理变量后，计算目标函数和约束条件以及
−6
−6
别为密度和杨氏模量插值函数的惩罚因子。 对应的敏度信息，最后通过移动渐进线法的全局求

解器（GCMMA）进行迭代求解 [36] 。GCMMA 求解中
ˆ ˆ ˆ
ξ ξ=Lξ ξ e ξ b ξ d
设 计 变 量 的 变 化 幅 值 不 超过 0.1， 初 始 投 影 强 度
β = 1，每经过 30 次外部迭代时，投影强度 β 变为原来
的 1.587 倍，直到 β 达到 64 时停止优化迭代过程，此
0: 空腔 固体: 1 时一共经历 300 代。

图 7 设计变量投影到物理变量的处理过程 图 8 中 右 侧 示 意 图 展 现 了 拓 扑 优 化 的 计 算 过
程，首先基于流固耦合传递矩阵法进行声学性能计
Fig. 7 Projection process of mapping design variables to
physical variables 算，然后基于拓扑优化算法在刚度、密度约束条件