基于密度法的多层约束阻尼结构拓扑优化设计

Topology Optimization Design of Multi‑layer Constrained Damping Structure Based on Density Method

摘要

针对高速织机等大幅面结构的抑振降噪问题，通过约束阻尼结构的本构关系，探索多阻尼层、约束层之间的数学模型；基于有限元方法建立多层约束阻尼板的动力学模型，利用拓扑优化理论，分析动力学模型中力学矩阵与单元密度设计变量间的耦合关系，建立多层约束阻尼板材料插值模型；基于密度法搭建多层约束阻尼板结构优化设计框架，并以约束阻尼材料的表面敷设总量为约束，以模态损耗因子的倒数最小化为优化目标，推导优化设计变量的灵敏度；最后以悬臂梁与两端固定梁为例对所提方法的有效性进行验证，结果表明：随着阻尼材料体积分数增加，模态损耗因子显著提高；在50%的体积约束下，两端固定梁的损耗因子从0.262提升至0.349，悬臂梁的损耗因子从0.138提升至0.215。

Abstract

To address the vibration suppression and noise reduction in large-scale structures such as high-speed looms, this study explores the mathematical model for multiple damping layers and constraining layers through the constitutive relationship of constrained damping structures. The dynamic model of the multi-layer constrained damping plate was then developed using the finite element method. By applying topology optimization theory, the coupling relationship between the mechanical matrix and the element density design variables in the dynamic model was analyzed, leading to the establishment of a material interpolation model for the multi-layer constrained damping plate.A structural optimization design framework for the multi-layer constrained damping plate structure was constructed based on the density method, with the optimization objective set to minimize the reciprocal of the modal loss factor, subject to a constraint on the total volume of the damping material. The sensitivity of the design variables was derived analytically.Finally, the effectiveness of the proposed method was verified using a cantilever beam and the fixed-fixed beam as examples. The results showed that the modal loss factor increased significantly with the increase of the damping material volume fraction. Under a 50% volume constraint, the loss factor of the fixed-fixed beam increased from 0.262 to 0.349, and that of the cantilever beam increased from 0.138 to 0.215.

Graphical abstract

关键词

减振 / 多层约束阻尼 / 材料插值模型

Key words

vibration reduction / multi-layer constrained damping / material interpolation model

随着纺织产业的快速发展，对纺织装备的抑振能力提出了更高的要求^［1］。尽管传统隔振技术能够实现基本的减振降噪功能，但其固有的局限性（如增加设备重量、恶化散热性能等）使其难以有效应对板壳结构振动这一复杂工程问题^［2］，而拓扑优化技术和阻尼材料的非线性特性在提高结构减振性能的同时还可以节省材料成本^［3］。因此，依托先进手段完成“安静地织布”已经成为纺织行业的重要研究方向^［4］。

在阻尼结构分析中，有限元方法展现了显著优势。Lu等^［5］采用有限元法对约束阻尼的夹层梁进行系统分析，证实了该方法在处理复杂阻尼结构中的适用性；Johnson等^［6］通过建立夹层梁有限元模型，在对比研究上运用了模态阻尼能量法，其数值解与实验结果的良好吻合验证了有限元公式的准确性。在阻尼敷设优化方面，Parthasarthy等^［7］借助有限元分析揭示了自由阻尼层局部敷设对矩形板模态参数的影响规律；Kumar等^［8］基于模态应变能法成功实现了阻尼片最佳敷设位置的确定；顾赛克等^［9］通过多参数组合研究，系统阐明了约束阻尼板结构的力学特性与减振机理，在多层约束阻尼结构的优化设计方面给予了重要参考。

在板壳减振性能上，阻尼材料的空间分布成为了关键影响因素^［10］。为此，学者们引入拓扑优化技术开展阻尼材料布局的优化研究^［11］。Kim等^［12］通过对比研究发现，以模态阻尼比最大化为目标的拓扑优化方法，获得的阻尼材料分布方案能产生更优异的模态损耗因子，显著优于传统设计方法；Qin等^［13］研究了多层阻尼材料层数对复合材料损耗因子的影响，发现约束阻尼结构的复合损耗因子主要受约束层的厚度、密度和模量的影响；Wan等^［14］通过系统实验研究了多层约束阻尼结构的减振特性，发现通过优化阻尼层的数量、厚度及空间排布等参数组合，可显著降低结构固有频率并有效提升损耗因子，证实了多层阻尼结构在振动控制方面的优越性能；Tobias等^［15］创新性地提出了2种阻尼结构拓扑优化策略——完全非线性方法和简化方法，成功实现了在考虑结构损伤的条件下，针对预定载荷工况的阻尼材料最优分布设计。

本研究基于多层复合结构的动力学特性分析，结合拓扑优化方法，构建了以模态损耗因子最大化为目标的多层板壳结构设计方法。对比单层板壳结构，本方法进一步提高了板壳结构的固有特性和吸能减振的能力，为纺织类板壳结构提供更好的选择方案。

1 多层约束阻尼板有限元建模

1.1 约束阻尼板

为了更加直观地反映多层约束阻尼板壳结构，本文通过两层约束阻尼的五层板来推导多层板的建模过程。为提升数值计算效率，在建立理论模型时暂未纳入频率、温度等可能影响材料特性的参数变量；同时，在优化过程中仅选取频率响应的实部作为计算依据。

两层约束阻尼结构示意图如图1所示。图1中u_{\mathrm{p}}、u_{{\mathrm{v}}_{\mathrm{1}}}、u_{{\mathrm{c}}_{\mathrm{1}}}、u_{{\mathrm{v}}_{\mathrm{2}}}、u_{{\mathrm{c}}_{\mathrm{2}}}分别表示基层、第一层阻尼层、第一层约束层、第二层阻尼层和第二层约束层沿x轴方向的形变；h\mathrm{p}、h{\mathrm{v}}_{\mathrm{1}}、h{\mathrm{c}}_{\mathrm{1}}、h{\mathrm{v}}_{\mathrm{2}}、h{\mathrm{c}}_{\mathrm{2}}分别表示基层、第一层阻尼层、第一层约束层、第二层阻尼层和第二层约束层的厚度；ɑ₁表示基层在x方向上的转角；ɑ₅表示第二层约束层在x方向上的转角；w为板的横向位移。

1.2 模型假设

相较于自由阻尼层，约束阻尼结构在动态载荷作用下呈现更复杂的变形行为。为研究方便，本文在建立多层约束阻尼板的动力学模型时，采用了薄板理论和线性假设，具体如下^［16］：

1） 几何假设：板的厚度远小于其最小平面尺寸；变形前垂直于中面的法线在变形后仍然保持为直线，且垂直于变形后的中面。

2） 变形假设：刚性的金属层采用Kirchhoff薄板理论（这是因为金属层厚度与其最小平面尺寸比\le \mathrm{1}/\mathrm{20}），并忽略其横向剪切变形；黏弹性阻尼层考虑其关键的面外剪切变形（这是因为阻尼层的耗能主要来自面外剪切变形），则采用Mindlin假设。

3） 应力-应变假设：垂直于板面的正应力{\sigma }_z远小于面内应力{\sigma }_x、{\sigma }_y，可忽略不计；面内应力{\sigma }_x、{\sigma }_y\mathrm{、}{\tau }_{xy}沿厚度线性分布，与弯曲力矩和扭矩相关。

4） 小变形与线性假设：板的横向位移w远小于板的厚度h，因此几何非线性可忽略；各层材料在目标工作频率内可近似视为线性、各向同性，并服从胡克定律。

1.3 各层形变的几何关系

基于以上模型假设，两层约束阻尼板沿着x轴方向的形变可以表达为：

（1）

式中：u_i代表第i层沿x方向的形变，其中i=\mathrm{p}\mathrm{、}{\mathrm{v}}_{\mathrm{1}}\mathrm{、}{\mathrm{c}}_{\mathrm{1}}\mathrm{、}{\mathrm{v}}_{\mathrm{2}}\mathrm{、}{\mathrm{c}}_{\mathrm{2}}分别表示基层、第一层阻尼层、第一层约束层、第二层阻尼层和第二层约束层；d表示不同材料层的表面与中性层的距离；u_{\mathrm{0}}表示各层中性层面沿x方向的形变。

则第一层阻尼层以及第二层阻尼层上、下表面的形变可表示为：

\left\{\begin{array}{l}{u}_{\mathrm{1}}=u_{\mathrm{p}}+\frac{h_{\mathrm{p}}}{\mathrm{2}}\frac{\partial w}{\partial x}\\ u_{\mathrm{2}}=u_{{\mathrm{c}}_{\mathrm{1}}}-\frac{h_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{2}}\frac{\partial w}{\partial x}\\ u_{\mathrm{3}}=u_{{\mathrm{c}}_{\mathrm{1}}}+\frac{h_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{2}}\frac{\partial w}{\partial x}\\ u_{\mathrm{4}}=u_{{\mathrm{c}}_{\mathrm{2}}}-\frac{h_{{\mathrm{c}}_{\mathrm{2}}}}{\mathrm{2}}\frac{\partial w}{\partial x}\end{array}\right.（2）

式中：u_{\mathrm{1}}、u_{\mathrm{2}}分别为第一层阻尼层下表面和上表面的形变；u_{\mathrm{3}}、u_{\mathrm{4}}分别为第二层阻尼层下表面和上表面的形变。

结合图1的单元位移形变关系，进一步得出两个阻尼层沿着x、y轴方向的形变为：

\left\{\begin{array}{l}{u}_{{\mathrm{v}}_{\mathrm{1}}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(u_{\mathrm{p}}+u_{{\mathrm{c}}_{\mathrm{1}}}\right)+\frac{h_{\mathrm{p}}-h_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{2}}\frac{\partial w}{\partial x}\right]\\ v_{{\mathrm{v}}_{\mathrm{1}}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(v_{\mathrm{p}}+v_{{\mathrm{c}}_{\mathrm{1}}}\right)+\frac{h_{\mathrm{p}}-h_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{2}}\frac{\partial w}{\partial y}\right]\end{array}\right.（3）

\left\{\begin{array}{l}{u}_{{\mathrm{v}}_{\mathrm{2}}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(u_{{\mathrm{c}}_{\mathrm{1}}}+u_{{\mathrm{c}}_{\mathrm{2}}}\right)+\frac{h_{{\mathrm{c}}_{\mathrm{1}}}-h_{{\mathrm{c}}_{\mathrm{2}}}}{\mathrm{2}}\frac{\partial w}{\partial x}\right]\\ v_{{\mathrm{v}}_{\mathrm{2}}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(v_{{\mathrm{c}}_{\mathrm{1}}}+v_{{\mathrm{c}}_{\mathrm{2}}}\right)+\frac{h_{{\mathrm{c}}_{\mathrm{1}}}-h_{{\mathrm{c}}_{\mathrm{2}}}}{\mathrm{2}}\frac{\partial w}{\partial y}\right]\end{array}\right.（4）

式中：v_{{\mathrm{v}}_{\mathrm{1}}}\mathrm{、}{v}_{{\mathrm{v}}_{\mathrm{2}}}分别为第一层阻尼层、第二层阻尼层沿y轴方向的形变。

1.4 形变模式和形函数

本研究中约束阻尼结构的有限元模型，采用四节点二维单元进行离散化处理。每个节点配置了9个自由度，分别为：1） 基层、第一层约束层、第二层约束层在x、y轴方向的形变u_{\mathrm{p}}、u_{{\mathrm{c}}_{\mathrm{1}}}、u_{{\mathrm{c}}_{\mathrm{2}}}\mathrm{、}{v}_{\mathrm{p}}、v_{{\mathrm{c}}_{\mathrm{1}}}、v_{{\mathrm{c}}_{\mathrm{2}}}； 2） 板的横向位移w； 3） 法线绕x、 y轴的转动{\theta }_x=\frac{\partial w}{\partial y}、{\theta }_y=-\frac{\partial w}{\partial x}。因此，单元节点第b个节点的形变向量可以表示为：

\left\{\begin{array}{l}{\mathit{q}}_b={\left[u_{\mathrm{p},b},v_{\mathrm{p},b},u_{{\mathrm{c}}_{\mathrm{1}},b},v_{{\mathrm{c}}_{\mathrm{1}},b},u_{{\mathrm{c}}_{\mathrm{2}},b},v_{{\mathrm{c}}_{\mathrm{2}},b},w,{\theta }_x,{\theta }_y\right]}^{\mathrm{T}}\\ b=\mathrm{1,2},\mathrm{3,4}\end{array}\right.（5）

式中： q_b 表示第b个节点的形变向量。

单元节点形变向量\mathit{q}可表示为：

\mathit{q}=\left[q_{\mathrm{1}},q_{\mathrm{2}},q_{\mathrm{3}},q_{\mathrm{4}}\right]（6）

式中：q₁、q₂、q₃、q₄分别表示约束阻尼结构的4个节点。

为定义单元内的形变场，采用带未知系数的多项式构造形变函数。为了确保形变模式在x、y方向上的对称性，插值函数设定为：

\left\{\begin{array}{l}{u}_{\mathrm{p}}={\alpha }_{\mathrm{1}}+{\alpha }_{\mathrm{2}}x+{\alpha }_{\mathrm{3}}y+{\alpha }_{\mathrm{4}}xy\\ v_{\mathrm{p}}={\alpha }_{\mathrm{5}}+{\alpha }_{\mathrm{6}}x+{\alpha }_{\mathrm{7}}y+{\alpha }_{\mathrm{8}}xy\\ u_{{\mathrm{c}}_{\mathrm{1}}}={\alpha }_{\mathrm{9}}+{\alpha }_{\mathrm{10}}x+{\alpha }_{\mathrm{11}}y+{\alpha }_{\mathrm{12}}xy\\ v_{{\mathrm{c}}_{\mathrm{1}}}={\alpha }_{\mathrm{13}}+{\alpha }_{\mathrm{14}}x+{\alpha }_{\mathrm{15}}y+{\alpha }_{\mathrm{16}}xy\\ u_{{\mathrm{c}}_{\mathrm{2}}}={\alpha }_{\mathrm{17}}+{\alpha }_{\mathrm{18}}x+{\alpha }_{\mathrm{19}}y+{\alpha }_{\mathrm{20}}xy\\ v_{{\mathrm{c}}_{\mathrm{2}}}={\alpha }_{\mathrm{21}}+{\alpha }_{\mathrm{22}}x+{\alpha }_{\mathrm{23}}y+{\alpha }_{\mathrm{24}}xy\\ \omega ={\alpha }_{\mathrm{25}}+{\alpha }_{\mathrm{26}}x+{\alpha }_{\mathrm{27}}y+{\alpha }_{\mathrm{28}}{x}^{\mathrm{2}}+{\alpha }_{\mathrm{29}}xy+\\ {\alpha }_{\mathrm{30}}{y}^{\mathrm{2}}+{\alpha }_{\mathrm{31}}{x}^{\mathrm{3}}+{\alpha }_{\mathrm{32}}{x}^{\mathrm{2}}y+{\alpha }_{\mathrm{33}}x{y}^{\mathrm{2}}+\\ {\alpha }_{\mathrm{34}}{y}^{\mathrm{3}}+{\alpha }_{\mathrm{35}}{x}^{\mathrm{3}}y+{\alpha }_{\mathrm{36}}x{y}^{\mathrm{3}}\\ {\theta }_x={\alpha }_{\mathrm{27}}+{\alpha }_{\mathrm{29}}x+\mathrm{2}{\alpha }_{\mathrm{30}}y+{\alpha }_{\mathrm{32}}{x}^{\mathrm{2}}+\\ \mathrm{2}{\alpha }_{\mathrm{33}}xy+\mathrm{3}{\alpha }_{\mathrm{34}}{y}^{\mathrm{2}}+{\alpha }_{\mathrm{35}}{x}^{\mathrm{3}}+\mathrm{3}{\alpha }_{\mathrm{36}}x{y}^{\mathrm{2}}\\ {\theta }_y={\alpha }_{\mathrm{26}}-\mathrm{2}{\alpha }_{\mathrm{28}}x-{\alpha }_{\mathrm{29}}y-\mathrm{3}{\alpha }_{\mathrm{31}}{x}^{\mathrm{2}}-\\ \mathrm{2}{\alpha }_{\mathrm{32}}xy-{\alpha }_{\mathrm{33}}{y}^{\mathrm{2}}-\mathrm{3}{\alpha }_{\mathrm{35}}{x}^{\mathrm{2}}y-{\alpha }_{\mathrm{36}}{y}^{\mathrm{3}}\end{array}\right.（7）

式中：{\alpha }_i,i=\mathrm{1},\cdots ,\mathrm{36}为待确定系数。

将式（7）简化成矩阵形式，可以表示为：

\complement \mathit{\alpha }={\mathit{q}}^{\mathrm{e}}=(u_{\mathrm{p}},v_{\mathrm{p}},u_{{\mathrm{c}}_{\mathrm{1}}},v_{{\mathrm{c}}_{\mathrm{1}}},u_{{\mathrm{c}}_{\mathrm{2}}},v_{{\mathrm{c}}_{\mathrm{2}}},w,{\theta }_x,{\theta }_y)（8）

式中：矩阵\complement的维度为36×36，其元素由单元节点的坐标值决定。

通过求逆运算可解出待定系数\mathit{\alpha }，即

\mathit{\alpha }={\complement }^{-\mathrm{1}}{\mathit{q}}^{\mathrm{e}}（9）

将解得的待定系数代回式（7），即可获得约束阻尼板内任意位置形变的插值表达式：

\begin{array}{l}{\left[u_{\mathrm{p}},v_{\mathrm{p}},u_{{\mathrm{c}}_{\mathrm{1}}},v_{{\mathrm{c}}_{\mathrm{1}}},u_{{\mathrm{c}}_{\mathrm{2}}},v_{{\mathrm{c}}_{\mathrm{2}}},w,{\theta }_x,{\theta }_y\right]}^{\mathrm{T}}=\\ \left[{\mathit{N}}_x^{\mathrm{p}},{\mathit{N}}_y^{\mathrm{p}},{\mathit{N}}_x^{{\mathrm{c}}_{\mathrm{1}}},{\mathit{N}}_y^{{\mathrm{c}}_{\mathrm{1}}},{\mathit{N}}_x^{{\mathrm{c}}_{\mathrm{2}}},{\mathit{N}}_y^{{\mathrm{c}}_{\mathrm{2}}},{\mathit{N}}_w,{\mathit{N}}_x^{\theta },{\mathit{N}}_y^{\theta }\right]{\mathit{q}}^{\mathrm{e}}\end{array}（10）

式中：{\mathit{N}}_x^{\mathrm{p}}\mathrm{、}{\mathit{N}}_y^{\mathrm{p}}\mathrm{、}{\mathit{N}}_x^{{\mathrm{c}}_{\mathrm{1}}}\mathrm{、}{\mathit{N}}_y^{{\mathrm{c}}_{\mathrm{1}}}\mathrm{、}{\mathit{N}}_x^{{\mathrm{c}}_{\mathrm{2}}}\mathrm{、}{\mathit{N}}_y^{{\mathrm{c}}_{\mathrm{2}}}\mathrm{、}{\mathit{N}}_w\mathrm{、}{\mathit{N}}_x^{\theta }\mathrm{、}{\mathit{N}}_y^{\theta }分别是u_{\mathrm{p}}\mathrm{、}{v}_{\mathrm{p}}\mathrm{、}{u}_{{\mathrm{c}}_{\mathrm{1}}}\mathrm{、}{v}_{{\mathrm{c}}_{\mathrm{1}}}\mathrm{、}{u}_{{\mathrm{c}}_{\mathrm{2}}}\mathrm{、}{v}_{{\mathrm{c}}_{\mathrm{2}}}\mathrm{、}w\mathrm{、}{\theta }_x\mathrm{、}{\theta }_y对应的形函数矩阵，每个矩阵的维度均为9×36。

式（10）对应于u_{{\mathrm{v}}_{\mathrm{1}}}\mathrm{、}{v}_{{\mathrm{v}}_{\mathrm{1}}}\mathrm{、}{u}_{{\mathrm{v}}_{\mathrm{2}}}\mathrm{、}{v}_{{\mathrm{v}}_{\mathrm{2}}}的形函数N_{u_{\mathrm{v}\mathrm{1}}}、N_{v_{\mathrm{v}\mathrm{1}}}、N_{u_{\mathrm{v}\mathrm{2}}}、N_{v_{\mathrm{v}\mathrm{2}}}可以表示为：

\left\{\begin{array}{l}{N}_{u_{\mathrm{v}\mathrm{1}}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(N_x^{\mathrm{p}}+N_x^{{\mathrm{c}}_{\mathrm{1}}}\right)-\frac{h_{\mathrm{p}}-h_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{2}}{N}_y^{\theta }\right]\\ N_{v_{\mathrm{v}\mathrm{1}}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(N_y^{\mathrm{p}}+N_y^{{\mathrm{c}}_{\mathrm{1}}}\right)+\frac{h_{\mathrm{p}}-h_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{2}}{N}_x^{\theta }\right]\\ N_{u_{\mathrm{v}\mathrm{2}}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(N_x^{{\mathrm{c}}_{\mathrm{1}}}+N_x^{{\mathrm{c}}_{\mathrm{2}}}\right)-\frac{h_{{\mathrm{c}}_{\mathrm{1}}}-h_{{\mathrm{c}}_{\mathrm{2}}}}{\mathrm{2}}{N}_y^{\theta }\right]\\ N_{v_{\mathrm{v}\mathrm{2}}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(N_y^{{\mathrm{c}}_{\mathrm{1}}}+N_y^{{\mathrm{c}}_{\mathrm{2}}}\right)+\frac{h_{{\mathrm{c}}_{\mathrm{1}}}-h_{{\mathrm{c}}_{\mathrm{2}}}}{\mathrm{2}}{N}_x^{\theta }\right]\end{array}\right.（11）

1.5 本构方程

基于式（1）~（4）的各层形变表达式与弹性力学理论，可推导出约束阻尼结构的单元应变场为：

\left\{\begin{array}{l}{\epsilon }_{lx}=\frac{\partial u}{\partial x}=\frac{\partial u_l}{\partial x}+z\frac{{\partial }^{\mathrm{2}}w}{\partial x^{\mathrm{2}}}\\ {\epsilon }_{ly}=\frac{\partial v}{\partial y}=\frac{\partial v_l}{\partial y}+z\frac{{\partial }^{\mathrm{2}}w}{\partial y^{\mathrm{2}}}\\ {\gamma }_{lxy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=\frac{\partial u_l}{\partial y}+\frac{\partial v_l}{\partial x}+\mathrm{2}z\frac{{\partial }^{\mathrm{2}}w}{\partial xy}\end{array}\right.（12）

式中：{\epsilon }_{lx}\mathrm{、}{\epsilon }_{ly}分别表示第l层沿x、y向的正应变，l=\mathrm{p},{\mathrm{v}}_{\mathrm{1}},{\mathrm{c}}_{\mathrm{1}},\cdots ,{\mathrm{v}}_n,{\mathrm{c}}_n为基层、第一层阻尼层、第一层约束层、……、第n层阻尼层、第n层约束层的几何特性（n表示阻尼层及其对应约束层的总层数）；{\gamma }_{lxy}表示第l层沿x‒y面的剪切应变；u_l\mathrm{、}{v}_l分别表示第l层沿x轴、y轴方向的面内形变。

各层对应的正应力可表示为：

\left\{\begin{array}{l}{\sigma }_{lx}=\frac{E_l}{\mathrm{1}-{\upsilon }_l^{\mathrm{2}}}\left({\epsilon }_{lx}+{\upsilon }_l{\epsilon }_{ly}\right)\\ {\sigma }_{ly}=\frac{E_l}{\mathrm{1}-{\upsilon }_l^{\mathrm{2}}}\left({\epsilon }_{ly}+{\upsilon }_l{\epsilon }_{lx}\right)\end{array}\right.（13）

式中：{\sigma }_{lx}\mathrm{、}{\sigma }_{ly}分别表示第l层沿x、y向的正应力；{\upsilon }_l表示第l层的泊松比。

各层对应的剪切应力可表示为：

{\tau }_{lxy}=G_l{\gamma }_{lxy}=\frac{E_l}{\mathrm{1}-{\upsilon }_l^{\mathrm{2}}}{\gamma }_{lxy}（14）

式中：{\tau }_{lxy}表示第l层沿x‒y面的剪切应力；E_l\mathrm{、}{G}_l分别表示第l层的弹性模量和剪切模量。

在后续能量分析中，基于板壳理论并忽略转动惯量的影响，采用能量法建立约束阻尼板的动能与应变能表达式。鉴于各层动能具有相同形式，因此可统一表达为：

\begin{array}{l}{E}_l^{\mathrm{e}}=\frac{\mathrm{1}}{\mathrm{2}}{\rho }_l{g}_l\iiint \left[{\left(\frac{\partial u_l}{\partial t}\right)}^{\mathrm{2}}+{\left(\frac{\partial v_l}{\partial t}\right)}^{\mathrm{2}}+{\left(\frac{\partial w_l}{\partial t}\right)}^{\mathrm{2}}\right]\mathrm{d}V=\\ \frac{\mathrm{1}}{\mathrm{2}}{\left({\mathit{q}}^{\mathrm{e}}\right)}^{\mathrm{T}}{\mathit{M}}_l^{\mathrm{e}}{\mathit{q}}^{\mathrm{e}}\end{array}（15）

式中：E_l^{\mathrm{e}}表示第l层的动能；{\rho }_l为第l层材料的密度；g_l 表示第l层材料的几何系数；w_l表示第l层的横向位移； M_l 表示第l层的单元质量矩阵；V为单元体积。

基于前节推导的不同材料层形函数，将式（7）、（10）代入式（15）可得：

\begin{array}{l}{E}_l^{\mathrm{e}}=\frac{\mathrm{1}}{\mathrm{2}}{\rho }_l{g}_l{\left({\mathit{q}}^{\mathrm{e}}\right)}^{\mathrm{T}}\times \\ {\int }_{-a}^a{\int }_{-b}^b\left(\begin{array}{l}{\left[{\mathit{N}}_{u_l}\right]}^{\mathrm{T}}\left[{\mathit{N}}_{u_l}\right]+\\ {\left[{\mathit{N}}_{v_l}\right]}^{\mathrm{T}}\left[{\mathit{N}}_{v_l}\right]+{\left[{\mathit{N}}_{w_l}\right]}^{\mathrm{T}}\left[{\mathit{N}}_{w_l}\right]\end{array}\right)\mathrm{d}x\mathrm{d}y{\mathit{q}}^{\mathrm{e}}\end{array}（16）

式中：a、b表示有限元单元在x、y轴方向的半长度；{\mathit{N}}_{u_l}\mathrm{、}{\mathit{N}}_{v_l}分别表示第l层沿x、y轴方向形变的形函数矩阵；{\mathit{N}}_{w_l}表示第l层横向位移w的形函数矩阵。

根据式（16），不同材料层的质量矩阵{\mathit{M}}_l^{\mathrm{e}}可表示为：

\begin{array}{l}{\mathit{M}}_l^{\mathrm{e}}={\rho }_l{g}_l\times \\ {\int }_{-a}^a{\int }_{-b}^b\left(\begin{array}{l}{\left[{\mathit{N}}_{u_l}\right]}^{\mathrm{T}}\left[{\mathit{N}}_{u_l}\right]+\\ {\left[{\mathit{N}}_{v_l}\right]}^{\mathrm{T}}\left[{\mathit{N}}_{v_l}\right]+{\left[{\mathit{N}}_{w_l}\right]}^{\mathrm{T}}\left[{\mathit{N}}_{w_l}\right]\end{array}\right)\mathrm{d}x\mathrm{d}y\end{array}（17）

进而，整体单元质量矩阵{\mathit{M}}_O^{\mathrm{e}}可表示为：

{\mathit{M}}_O^{\mathrm{e}}={\displaystyle \sum _{n=\mathrm{1}}^O}\left({\mathit{M}}_p^e+{\mathit{M}}_{v_n}^e+{\mathit{M}}_{c_n}^e\right)（18）

式中：{\mathit{M}}_{\mathrm{p}}^{\mathrm{e}}为基层质量矩阵；{\mathit{M}}_{{\mathrm{v}}_n}^{\mathrm{e}}表示第n层阻尼层的质量矩阵；{\mathit{M}}_{{\mathrm{c}}_n}^{\mathrm{e}}表示第n层约束层的质量矩阵；O表示约束层与阻尼层的总套数。

单元在不同材料层上的势能T_l^{\mathrm{e}}可表示为：

\begin{array}{l}{T}_l^{\mathrm{e}}=\frac{\mathrm{1}}{\mathrm{2}}\iiint \left({\sigma }_{lx}{\epsilon }_{lx}+{\sigma }_{ly}{\epsilon }_{ly}+{\tau }_{lxy}{\gamma }_{lxy}\right)\mathrm{d}V=\\ \frac{\mathrm{1}}{\mathrm{2}}{\left({\mathit{q}}^{\mathrm{e}}\right)}^{\mathrm{T}}{\mathit{K}}_l^{\mathrm{e}}{\mathit{q}}^{\mathrm{e}}\end{array}（19）

式中：{\mathit{K}}_l^{\mathrm{e}}表示不同材料层的单元刚度矩阵。

结合胡克定律及单元节点形变对应的形函数，将式（12）~（14）代入式（19），势能表达式转化为：

T_l^{\mathrm{e}}=\frac{\mathrm{1}}{\mathrm{2}}{\left({\mathit{q}}^{\mathrm{e}}\right)}^{\mathrm{T}}\left(h_l{\int }_{-a}^a{\int }_{-b}^b{{\mathit{B}}_l}^{\mathrm{T}}{\mathit{D}}_l{\mathit{B}}_l\mathrm{d}x\mathrm{d}y+\frac{h_l^{\mathrm{3}}}{\mathrm{12}}{\int }_{-a}^a{\int }_{-b}^b{{\mathit{B}}_l}^{\mathrm{T}}{\mathit{D}}_l{\mathit{B}}_l\mathrm{d}x\mathrm{d}y\right){\mathit{q}}^{\mathrm{e}}（20）

式中：h_l表示第l层的厚度；{\mathit{B}}_l表示第l层应变转换矩阵；{\mathit{D}}_l表示第l层弹性常数矩阵。

由此导出不同材料层单元刚度矩阵{\mathit{K}}_l^{\mathrm{e}}为：

{\mathit{K}}_l^{\mathrm{e}}=h_l{\int }_{-a}^a{\int }_{-b}^b{{\mathit{B}}_l}^{\mathrm{T}}{\mathit{D}}_l{\mathit{B}}_l\mathrm{d}x\mathrm{d}y+\frac{h_l^{\mathrm{3}}}{\mathrm{12}}{\int }_{-a}^a{\int }_{-b}^b{{\mathit{B}}_l}^{\mathrm{T}}{\mathit{D}}_l{\mathit{B}}_l\mathrm{d}x\mathrm{d}y（21）

应变转换矩阵{\mathit{B}}_l、弯曲应变转换矩阵 B 和弹性常数矩阵{\mathit{D}}_l分别表征各层的几何材料特性，具体定义如下：

\left\{\begin{array}{l}{\mathit{B}}_l=\left[\begin{array}{cc}\frac{\partial N_{u_l}}{\partial x}& \mathrm{0}\\ \mathrm{0}& \frac{\partial N_{v_l}}{\partial y}\\ \frac{\partial N_{u_l}}{\partial y}& \frac{\partial N_{v_l}}{\partial x}\end{array}\right]\\ \mathit{B}=\left[\begin{array}{ccc}\frac{{\partial }^{\mathrm{2}}{N}_w}{\partial x^{\mathrm{2}}}& \frac{{\partial }^{\mathrm{2}}{N}_w}{\partial y^{\mathrm{2}}}& \frac{{\partial }^{\mathrm{2}}{N}_w}{\partial xy}\end{array}\right]\\ {\mathit{D}}_l=\frac{E_l}{\mathrm{1}-{{\mu }_l}^{\mathrm{2}}}\left(\begin{array}{ccc}\mathrm{1}& {\mu }_l& \mathrm{0}\\ {\mu }_l& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}-{\mu }_l}{\mathrm{2}}\end{array}\right)\end{array}\right.（22）

通过将约束阻尼层视为线性弹性材料，第n层阻尼层剪切应变产生的弹性势能T_{{\mathrm{v}}_n}^{\mathrm{e}}可以表示为：

T_{{\mathrm{v}}_n}^{\mathrm{e}}=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \underset{V}{\overset{}{\iiint }}}\left(G_{{\mathrm{v}}_n}{\beta }_{x_n}^{\mathrm{2}}+G_{{\mathrm{v}}_n}{\beta }_{y_n}^{\mathrm{2}}\right)\mathrm{d}V=\frac{\mathrm{1}}{\mathrm{2}}{\left({\mathit{q}}^{\mathrm{e}}\right)}^{\mathrm{T}}{\mathit{K}}_{{\beta }_{\mathrm{v}n}}^{\mathrm{e}}{\mathit{q}}^{\mathrm{e}}（23）

式中；G_{{\mathrm{v}}_n}为第n层阻尼层的剪切模量；{\beta }_{x_n}为第n层阻尼层沿x方向的剪切应变；{\beta }_{y_n}为第n层阻尼层沿y方向的剪切应变；{\mathit{K}}_{{\beta }_{\mathrm{v}n}}^{\mathrm{e}}表示第n层阻尼层的剪切刚度矩阵。

依据式（21）及该阻尼层对应的形函数式（11），构造第n层阻尼层的剪切刚度矩阵{\mathit{K}}_{{\beta }_{\mathrm{v}n}}^{\mathrm{e}}为：

{\mathit{K}}_{{\beta }_{\mathrm{v}n}}^{\mathrm{e}}=G_{{\mathrm{v}}_n}^{*}{h}_{{\mathrm{v}}_n}{\int }_{-a}^a{\int }_{-b}^b({\left[{\mathit{N}}_{{\beta }_{xn}}\right]}^{\mathrm{T}}{\mathit{N}}_{{\beta }_{xn}}+{\left[{\mathit{N}}_{{\beta }_{yn}}\right]}^{\mathrm{T}}{\mathit{N}}_{{\beta }_{yn}})\mathrm{d}x\mathrm{d}y（24）

式中：G_{{\mathrm{v}}_n}^{*}为第n层阻尼层的复剪切模量；h_{{\mathrm{v}}_n}表示第n层阻尼层的厚度；{\mathit{N}}_{{\beta }_{xn}}\mathrm{、}{\mathit{N}}_{{\beta }_{yn}}分别表示第n层阻尼层沿x、y方向剪切应变的形函数矩阵。

因此，整体单元刚度矩阵{\mathit{K}}_O^{\mathrm{e}}可综合表示为：

{\mathit{K}}_O^{\mathrm{e}}={\displaystyle \sum _{n=\mathrm{1}}^O}\left({\mathit{K}}_p^e+{\mathit{K}}_{v_n}^e+{\mathit{K}}_{c_n}^e+{\mathit{K}}_{{\beta }_{vn}}^e\right)（25）

式中：{\mathit{K}}_{\mathrm{p}}^{\mathrm{e}}为基层刚度矩阵；{\mathit{K}}_{{\mathrm{v}}_n}^{\mathrm{e}}表示第n层阻尼层的刚度矩阵；{\mathit{K}}_{{\mathrm{c}}_n}^{\mathrm{e}}表示第n层约束层的刚度矩阵。

通过组装单元质量矩阵、单元刚度矩阵，最终得到整体结构的质量矩阵 M 、刚度矩阵 K 为

\left\{\begin{array}{l}\mathit{M}={\displaystyle \sum _{j=\mathrm{1}}^N}\left({\displaystyle \sum _{n=\mathrm{1}}^O}\left({\mathit{M}}_{\mathrm{p}}^{\mathrm{e}}+{\mathit{M}}_{{\mathrm{v}}_n}^{\mathrm{e}}+{\mathit{M}}_{{\mathrm{c}}_n}^{\mathrm{e}}\right)\right)\\ \mathit{K}={\displaystyle \sum _{j=\mathrm{1}}^N}\left({\displaystyle \sum _{n=\mathrm{1}}^O}\left({\mathit{K}}_{\mathrm{p}}^{\mathrm{e}}+{\mathit{K}}_{{\mathrm{v}}_n}^{\mathrm{e}}+{\mathit{K}}_{{\mathrm{c}}_n}^{\mathrm{e}}+{\mathit{K}}_{{\beta }_{\mathrm{v}n}}^{\mathrm{e}}\right)\right)\end{array}\right.（26）

式中：j为单元个数；N表示总阻尼单元个数。

1.6 模态应变能计算

黏弹性阻尼材料的应力-应变存在相位滞后效应，该特性可采用复模量模型进行表征。材料本构关系G^{*}可表述为：

G^{*}=G^{\text{'}}+\mathrm{i}{G}^{\text{'}\text{'}}=G^{\text{'}}\left(\mathrm{1}+\mathrm{i}\eta \right)（27）

式中：G^{\mathrm{\text{'}}}为阻尼材料的储能模量；G^{"}为阻尼本构关系的损耗模量；\eta =G^{\mathrm{\text{'}}\mathrm{\text{'}}}/G^{\mathrm{\text{'}}}为阻尼损耗因子；\mathrm{i}为虚数单位。

由此约束阻尼板的杨氏模量E^{*}可以定义为：

E^{*}=\mathrm{2}\left(\mathrm{1}+\mu \right)G^{*}（28）

式中：\mu为多层约束阻尼结构中各层材料的泊松比。

鉴于黏弹性材料的阻尼系数较小，可假设其阻尼系统与非阻尼系统的共振频率及振型基本吻合。因此，多层约束阻尼结构的动力学方程可建立为：

\mathit{M}\ddot{U}+\mathit{K}U=\mathrm{0}（29）

式中：U、\ddot{U}分别表示多层约束阻尼结构整体的位移和加速度；

式（29）解的形式可表示为u=\lambda \mathrm{s}\mathrm{i}\mathrm{n}\omega (t-t_{\mathrm{0}})（\lambda表示自由振动的特征值）。因此，黏弹性材料的模态方程为：

\left(\mathit{K}-\lambda \mathit{M}\right)U=\mathrm{0}（30）

由于两层阻尼层的材料属性都相同，则系统总模态应变能与第r阶模态阻尼比可表示为：

\left\{\begin{array}{l}{E}_{\mathrm{a}r}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\varphi }}^{\mathrm{T}}K\mathit{\varphi }\\ {\epsilon }_r=\frac{n_{{\mathrm{v}}_{\mathrm{1}}}{E}_{{\mathrm{v}}_{r\mathrm{1}}}+n_{{\mathrm{v}}_{\mathrm{2}}}{E}_{{\mathrm{v}}_{r\mathrm{2}}}}{E_{{\mathrm{p}}_r}+E_{{\mathrm{v}}_{r\mathrm{1}}}+E_{{\mathrm{v}}_{r\mathrm{2}}}+E_{{\mathrm{c}}_{r\mathrm{1}}}+E_{{\mathrm{c}}_{r\mathrm{1}}}}=\mathrm{2}{n}_{\mathrm{v}}\frac{{{\mathit{\varphi }}_r}^{\mathrm{T}}{\mathit{K}}_{\mathrm{v}}{\mathit{\varphi }}_r}{{{\mathit{\varphi }}_r}^{\mathrm{T}}\mathit{K}{\mathit{\varphi }}_r}\end{array}\right.（31）

式中：E_{\mathrm{a}r}表示第r阶模态的总应变能；\mathit{\varphi }为结构的特征向量，即模态振型；n_{\mathrm{v}}表示材料的损耗因子；E_{{\mathrm{p}}_r}\mathrm{、}{E}_{{\mathrm{v}}_{r\mathrm{1}}}\mathrm{、}{E}_{{\mathrm{v}}_{r\mathrm{2}}}\mathrm{、}{E}_{{\mathrm{c}}_{r\mathrm{1}}}\mathrm{、}{E}_{{\mathrm{c}}_{r\mathrm{1}}}分别为基层、第一层阻尼层、第二层阻尼层、第一层约束层、第二层约束层第r阶模态的应变能；{\mathit{\varphi }}_r为第r阶模态振型；{\mathit{K}}_{\mathrm{v}}为总的阻尼层模态应变能。

2 拓扑优化模型及灵敏度分析

2.1 拓扑优化模型

在多层约束阻尼结构的优化设计中，材料在结构中的空间分布状态直接影响其振动能量的耗散能力。传统的尺寸优化或形状优化方法仅能在预先设定的拓扑构型下调整参数，难以充分探索最优的材料布局方案。拓扑优化通过在设计域内对材料分布进行自由优化，可同时确定结构的几何形状与材料的排布方式，因而被广泛运用于轻量化与高性能结构设计中。

在拓扑优化方法中，惩罚性各向同性材料模型（solid isotropic material with penalization，SIMP）方法是一种成熟且计算效率较高的连续型材料插值策略，其核心思想可概括为：将每个有限单元的材料状态用一个连续设计变量x_j\in \left[x_{\mathrm{m}\mathrm{i}\mathrm{n}},\mathrm{1}\right]表示（x_j=\mathrm{1}时表示该单元完全为阻尼材料，x_j=x_{\mathrm{m}\mathrm{i}\mathrm{n}}时表示该单元几乎没有阻尼材料），并通过引入指数型惩罚因子将中间密度值推向0或1，从而获得接近离散型的材料分布；最后将结构的刚度矩阵、阻尼矩阵等力学参数通过插值形式表达，以保证材料属性与设计变量之间的数学连续性。

基于该方法，本研究以约束阻尼材料的表面敷设总量为约束，以模态损耗因子的倒数最小化为优化目标，实现结构减振吸能的设计目的。基于SIMP法的多层约束阻尼拓扑优化模型如下：

（32）

式中：x_j表示第j单元约束阻尼材料的相对密度；\Psi表示优化目标函数；{\xi }_r表示第r阶模态的损耗因子；V_j表示第j单元阻尼层的体积；x_{\mathrm{m}\mathrm{i}\mathrm{n}}为最小相对密度，本研究中设定x_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{0.001}。

2.2 灵敏度分析

拓扑优化过程中，需要计算目标函数对设计变量的偏导数进行灵敏度分析，以评估设计参数对目标函数梯度的影响，并基于该分析结果改进优化算法，以加速收敛过程并提升迭代效率。

结合多层阻尼材料的本构关系与优化模型，对阻尼层及约束层施加SIMP惩罚因子，可以推导出更新后的整体刚度矩阵\mathit{K}\mathit{\text{'}}与质量矩阵\mathit{M}\mathit{\text{'}}为：

\left\{\begin{array}{l}\mathit{K}\mathit{\text{'}}={\displaystyle \sum _{j=\mathrm{1}}^N}\left({\mathit{K}}_p^e+x_j^g\left({\displaystyle \sum _{n=\mathrm{1}}^O}\left({\mathit{K}}_{{\mathrm{v}}_n}^{\mathrm{e}}+{\mathit{K}}_{{\mathrm{c}}_n}^{\mathrm{e}}+{\mathit{K}}_{{\beta }_{\mathrm{v}n}}^{\mathrm{e}}\right)\right)\right)\\ \mathit{M}\mathit{\text{'}}={\displaystyle \sum _{j=\mathrm{1}}^N}\left({\mathit{M}}_p^e+x_j^q\left({\displaystyle \sum _{n=\mathrm{1}}^O}\left({\mathit{M}}_{{\mathrm{v}}_n}^{\mathrm{e}}+{\mathit{M}}_{{\mathrm{c}}_n}^{\mathrm{e}}\right)\right)\right)\end{array}\right.（33）

式中：g为刚度矩阵惩罚因子；q为质量矩阵惩罚因子；x_j^g、x_j^q分别表示第j个单元相对密度的g、q次幂。

将式（33）代入式（31），得到更新后的模态损耗因子{\xi }_r^{\text{'}}为

{\xi }_r^{\text{'}}=\frac{n_{{\mathrm{v}}_{\mathrm{1}}}{E}_{{\mathrm{v}}_{r\mathrm{1}}}+n_{{\mathrm{v}}_{\mathrm{2}}}{E}_{{\mathrm{v}}_{r\mathrm{2}}}}{E_{{\mathrm{p}}_r}+E_{{\mathrm{v}}_{r\mathrm{1}}}+E_{{\mathrm{v}}_{r\mathrm{2}}}+E_{{\mathrm{c}}_{r\mathrm{1}}}+E_{{\mathrm{c}}_{r\mathrm{1}}}}=\mathrm{2}{n}_{\mathrm{v}}\frac{{{\mathit{\varphi }}_r}^{\mathrm{T}}{\mathit{K}}_{\mathrm{v}}{\mathit{\varphi }}_r}{{{\mathit{\varphi }}_r}^{\mathrm{T}}\mathit{K}{\mathit{\varphi }}_r}（34）

经SIMP插值处理的整体刚度矩阵K\text{'}组装完成后，计算目标函数对设计变量的偏导数为：

\frac{\partial \Psi }{\partial x_j}=\frac{\left(\begin{array}{l}{{\mathit{\varphi }}_r}^{\mathrm{T}}{\mathit{K}}_{\mathrm{v}}{\mathit{\varphi }}_r\frac{\partial }{\partial x_j}\left({{\mathit{\varphi }}_r}^{\mathrm{T}}\mathit{K}{\mathit{\varphi }}_r\right)-\\ {{\mathit{\varphi }}_r}^{\mathrm{T}}\mathit{K}{\mathit{\varphi }}_r\frac{\partial }{\partial x_j}\left({{\mathit{\varphi }}_r}^{\mathrm{T}}{\mathit{K}}_{\mathrm{v}}{\mathit{\varphi }}_r\right)\end{array}\right)}{\mathrm{2}{\eta }_{\mathrm{v}}\cdot {\left({{\mathit{\varphi }}_r}^{\mathrm{T}}{\mathit{K}}_{\mathrm{v}}{\mathit{\varphi }}_r\right)}^{\mathrm{2}}}（35）

假设第r阶模态振型{\varphi }_r对x_j的导数忽略不计，则

\left\{\begin{array}{l}\frac{\partial }{\partial x_j}\left({{\mathit{\varphi }}_r}^{\mathrm{T}}\mathit{K}{\mathit{\varphi }}_r\right)={{\mathit{\varphi }}_r}^{\mathrm{T}}\frac{\partial \mathit{K}}{\partial x_j}{\mathit{\varphi }}_r\\ \frac{\partial }{\partial x_j}\left({{\mathit{\varphi }}_r}^{\mathrm{T}}{\mathit{K}}_{\mathrm{v}}{\mathit{\varphi }}_r\right)={{\mathit{\varphi }}_r}^{\mathrm{T}}\frac{\partial {\mathit{K}}_{\mathrm{v}}}{\partial x_j}{\varphi }_r\end{array}\right.（36）

将式（35）、（36）合并，得到最终结果为：

\frac{\partial \Psi }{\partial x_j}=\frac{\left(\begin{array}{l}\left({{\mathit{\varphi }}_r}^{\mathrm{T}}\frac{\partial \mathit{K}}{\partial x_j}{\mathit{\varphi }}_r\right)\left({{\mathit{\varphi }}_r}^{\mathrm{T}}{\mathit{K}}_{\mathrm{v}}{\mathit{\varphi }}_r\right)-\\ \left({{\mathit{\varphi }}_r}^{\mathrm{T}}\mathit{K}{\mathit{\varphi }}_r\right)\left({{\mathit{\varphi }}_r}^{\mathrm{T}}\frac{\partial {\mathit{K}}_{\mathrm{v}}}{\partial x_j}{\mathit{\varphi }}_r\right)\end{array}\right)}{\mathrm{2}{\eta }_{\mathrm{v}}\cdot {\left({{\mathit{\varphi }}_r}^{\mathrm{T}}{\mathit{K}}_{\mathrm{v}}{\mathit{\varphi }}_r\right)}^{\mathrm{2}}}（37）

对阻尼层及约束层施加SIMP惩罚因子后，经常存在数值不稳定的现象，需要采用Sigmund提出的灵敏度过滤技术^［17］，即通过设定最小过滤半径r_{\mathrm{m}\mathrm{i}\mathrm{n}}对目标单元实施邻域加权平均处理，具体方法为：1） 提升中心单元邻近区域的灵敏度权重；2） 抑制远离中心单元的低灵敏度影响；3） 原中心单元灵敏度值由过滤域加权均值替代。更新后的灵敏度表达式为：

\frac{\partial \Psi }{\partial x_j}=\frac{\mathrm{1}}{x_j{\displaystyle \sum _{f=\mathrm{1}}^{N\text{'}}}{\widehat{H}}_f}{\displaystyle \sum _{f=\mathrm{1}}^{N\text{'}}}({\widehat{H}}_f{x}_f\frac{\partial \Psi }{\partial x_f})（38）

其中，{\widehat{H}}_f=r_{\mathrm{m}\mathrm{i}\mathrm{n}}-∆d\left(m,z\right),z\in N\text{'},∆d\left(m,z\right)\le r_{\mathrm{m}\mathrm{i}\mathrm{n}}

式中：{\widehat{H}}_f为过滤函数，表示卷积算子；x_f表示过滤区域内第f单元的相对密度；N\text{'}表示过滤区域内的单元总数；∆d\left(m,z\right)为单元m、z的中心点距离。

为了让拓扑优化迭代稳定收敛，最终获得阻尼材料用量约束下模态损耗因子最优的多层约束阻尼结构构型，需要采用OC法迭代调整多层约束阻尼结构的设计参数。迭代公式为

x_{\mathrm{n}\mathrm{e}\mathrm{w}}=\left\{\begin{array}{l}\mathrm{m}\mathrm{a}\mathrm{x}\left(x_{\mathrm{m}\mathrm{i}\mathrm{n}},x-\epsilon \right),x{D}^s\le \mathrm{m}\mathrm{a}\mathrm{x}\left(x_{\mathrm{m}\mathrm{i}\mathrm{n}},x-\epsilon \right)\\ x{D}^s,\mathrm{m}\mathrm{a}\mathrm{x}\left(x_{\mathrm{m}\mathrm{i}\mathrm{n}},x-\epsilon \right)x{D}^s\le \mathrm{m}\mathrm{i}\mathrm{n}(\mathrm{1},x+\epsilon )\\ \mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{1},x+\epsilon \right),\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{1},x+\epsilon \right)\le x{D}^s\end{array}\right.（39）

其中 ， D=-\frac{\frac{\partial \Psi }{\partial x}}{K\frac{\partial V}{\partial x}}

式中：x、x_{\mathrm{n}\mathrm{e}\mathrm{w}}分别表示迭代更新前、后第j个单元的相对密度；\epsilon为移动极限，取值区间为（0，1），典型取值为0.3；s为迭代步数，记录优化过程循环次数；D为目标函数、体积约束对设计变量导数的比值；D^s表示优化准则系数D迭代步数的s次幂。

公式（39）迭代收敛判断的定义为：当连续两次迭代的相对密度差值的绝对值小于0.001时，即满足收敛条件。

3 优化算例

3.1 悬臂梁优化

为了深入研究多层约束阻尼结构的振动特性，对比不同体积分数下固有频率与模态损耗因子的变化，本节以图2所示的一端固定的悬臂梁多层约束阻尼结构为例进行优化设计。图2中多层约束阻尼结构基板材料采用铝板，密度为\mathrm{2}\mathrm{700}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，厚度为1.8 mm，弹性模量为\mathrm{6.89}\times {\mathrm{10}}^{\mathrm{4}} MPa，泊松比为0.3。在不考虑温度影响的情况下，两层阻尼层均采用恒温约束阻尼胶片，剪切模量为13.38 MPa，泊松比为0.495，密度为\mathrm{1}\mathrm{450}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，厚度为1.1 mm；两层约束层材料均采用铝板，其材料参数与基板完全相同，厚度为0.1 mm，过滤半径设为1.3。

图3为采用上述方法对一端固定的悬臂梁多层约束阻尼结构进行优化的结果。由图3可知，优化后的固有频率及模态损耗因子均随体积约束增大呈上升趋势，且均较优化前（体积约束为0）降低。虽然优化前的减振降噪效果较好，但考虑材料成本、设备重量、散热效果等因素，体积约束为50%的综合效果最好，即在维持结构固有特性的前提下，通过拓扑优化技术可同步实现阻尼材料减量与振动抑制。

为了分析不同网格划分对拓扑构型的影响，将多层约束阻尼拓扑的惩罚因子g设为3、q设为1，设计域网格分别划分为18×12、24×12、36×18。图4为上述状况下一阶模态损耗因子体积约束为50%时不同网格划分的拓扑优化构型。由图4可知，在整体优化构型图中虽然存在少量的中间密度单元，但材料边界清晰可辩；结构几何特征显著锐化，模态损耗因子呈现小幅提升；所有工况均保持快速收敛特性，但网格加密导致迭代次数增加，计算时长相应延长。

以图4b为例，对多层约束阻尼拓扑进行优化迭代，结果如图5所示。由图5可知，迭代曲线收敛平稳，模态损耗因子随迭代步数提高持续提升（由初始值迭代次数为2时的0.138增至迭代次数为50时的0.215）；中间密度单元占比随迭代推进显著降低，在第50步获得最优解，最终拓扑构型中中间密度单元趋近于0。该演化过程验证了算法的有效性与鲁棒性。

3.2 两端固定梁的优化

本节以图6所示的两端固定的多层约束阻尼结构为例，以一阶模态损耗因子最大化为目标进行优化设计。优化时，以多层约束阻尼使用量（体积）的40%、50%、60%为约束条件，刚度惩罚因子g和质量惩罚因子q分别设置为3和1，网格划分为48×12，得到不同体积约束下的优化图如图7所示。从图7可以看出，随着体积约束占比的提升，优化后的一阶模态损耗因子逐步提高；在体积约束50%、60%下，优化后的一阶模态损耗因子呈现非常清晰的拓扑构型，其中黑色部分为敷设的阻尼材料。

从优化效果与计算效率来看，60%体积约束下虽然模态损耗因子略高于50%的工况，但优化计算时长显著增加；50%体积约束下，结构已具备优异的减振吸能性能，且能有效控制计算成本，无需过度追求材料用量的提升，因此综合来看50%体积约束的优化方案效果最优。

图8对比了在体积分数为50%的最优工况下，单层与多层约束阻尼拓扑构型的演化。由图8可知，整个迭代过程收敛稳定；随优化步数增加，中间密度单元持续消减并渐近最优构型；相较于单层约束阻尼结构，多层约束阻尼结构损耗因子显著提升，吸能减振效果大大增强。

4 结论

基于多层约束阻尼结构的几何参数与材料参数，结合拓扑优化技术，建立以变密度法为核心的多层板壳性能优化框架；构建拓扑优化数学模型，提出灵敏度过滤计算策略；以模态损耗因子为目标函数，求解不同体积分数约束下的最优构型。以悬臂梁及两端固定梁为例进行试验研究，结果表明该方法能够显著提升结构的阻尼性能。

上述研究均在室温条件下进行，而实际装备运行中，板壳结构因高速运转必然产生温升；温度升高将显著改变阻尼层黏弹性材料的动态力学性能，进而削弱结构的吸能减振效果。因此，温度梯度下多层板壳的性能演化规律及优化设计方法将是下一步的研究重点。

参考文献

[1] 邱海飞. 综框的技术进展及研究方向［J］. 棉纺织技术， 2016，44（5）： 80-84.

[2] 刘棣华. 结构阻尼减振系统设计［J］. 强度与环境， 1983， 10（4）： 31-40， 55.

[3] 赵全鹏， 杨建成， 刘艳哲， 等. 经纱张力调节机构的拓扑优化设计与疲劳分析［J］. 现代纺织技术， 2023， 31（3）： 27-35.

[4] 邱海飞， 万宏强. 基于材料特性的综框专用件动态设计［J］. 丝绸， 2019， 56（7）： 59-64.

[5] LU Y P， EVERSTINE G C. More on finite element modeling of damped composite systems［J］. Journal of Sound and Vibration， 1980， 69（2）： 199-205.

[6] JOHNSON C D， KIENHOLZ D A. Finite element prediction of damping in structures with constrained viscoelastic layers［J］. Aiaa Journal， 1982， 20（9）： 1284-1290.DOI：10.2514/3.51190 .

[7] PARTHASARTHY G， GANESAN N， REDDY C V R. Study of unconstrained layer damping treatments applied to rectangular plates having central cutouts［J］. Computers & Structures， 1986， 23（3）： 433-443.

[8] KUMAR N， SINGH S P. Experimental study on vibration and damping of curved panel treated with constrained viscoelastic layer［J］. Composite Structures， 2010， 92（2）： 233-243.

[9] 顾赛克， 邓琼， 刘悦. 约束阻尼结构力学性能研究与参数优化［J］. 航空工程进展， 2021， 12（4）： 68-79， 89.

[10] 任晋宇， 徐静， 田波， 等. 阻尼覆盖面积对船舶约束阻尼板振动特性影响分析［J］. 噪声与振动控制， 2021， 41（1）： 122-126.

[11] 王兴锋， 张氢， 秦仙蓉， 等. 基于半定规划的结构轻量化设计［J］. 机械强度， 2022， 44（6）： 1365-1370.

[12] KIM S Y， MECHEFSKE C K， KIM I Y. Optimal damping layout in a shell structure using topology optimization［J］. Journal of Sound and Vibration， 2013， 332（12）： 2873-2883.

[13] QIN W Q， ZHANG J， LI J， et al. An inquiry into the uncertainty mechanism： the influence of the number of layers in multilayered damping materials on composite loss factor［J］. Polymer Composites， 2024， 45（10）： 9169-9180.

[14] WAN H C， LI Y N， ZHENG L. Vibration and damping analysis of a multilayered composite plate with a viscoelastic midlayer［J］. Shock and Vibration， 2016， 2016： 6354915.

[15] TOBIAS B， EMAD S， OLE S， et al.Topology optimization of damage‐resistant structures with a predefined load‐bearing capacity［J］.International Journal for Numerical Methods in Engineering，2021，123（4）：1114-1145.

[16] 陈姗姗. 基于变密度法的约束阻尼结构拓扑优化研究［D］. 桂林： 桂林电子科技大学， 2020.

[17] 张国锋， 徐雷， 王鑫， 等. 面向连续体结构拓扑优化的分区密度修正敏度过滤方法研究［J］.机械科学与技术， 2022， 41（11）： 1641-1649.