约束阻尼结构（constrained layer damping，CLD）主要依靠粘弹性阻尼材料的剪切变形消耗振动能量，实现主体板壳结构的振动噪声控制，具有结构简单、易于工程实现、减振降噪效果良好以及可靠性高等优点^[1]，广泛应用于舰船、航天飞行器、汽车等的承载板壳减振降噪设计中。但将粘弹性阻尼材料敷设于整个结构表面，不仅难以获得最优的减振降噪性能，而且会增加附加质量。因此，在轻量化设计背景下，为充分有效地利用材料，对约束阻尼结构进行优化设计，具有重要的工程意义。

约束阻尼结构的参数优化设计已有非常广泛的研究^[2-4]，多以模态损耗因子和结构的峰值响应为目标函数，对约束层和阻尼层材料的几何参数和材料参数等进行优化设计。参数优化多以表面全部敷设约束阻尼材料的结构为对象，这显然难以有效利用材料。随着连续结构拓扑优化的发展，拓扑优化方法被引入到约束阻尼结构的优化问题中。王明旭等^[5]借助固体各向同性微结构惩罚（solidisotropic microstructures with penalization，SIMP）模型，以模态阻尼比最大化为目标，对约束阻尼板壳结构进行拓扑寻优。Fang等^[6]以最小化谐振响应为目标，采用渐进优化方法（(evolutionary structural optimization，ESO）寻找板结构表面的约束阻尼材料的最优布局。窦松然等^[7]基于材料属性的有理近似（rational approximation of material prop-erties，RAMP）插值模型密度法，以模态损耗因子最大化为目标，对约束阻尼结构进行拓扑优化分析。Araújo等^[8]采用直接多元搜索法，以限定频率范围的模态损耗因子为目标，对约束阻尼夹层板进行多目标拓扑优化。Zhang等^[9]则以辐射声功率为目标函数，基于双向渐进优化（bi-direction evolutionary structural optimization，BESO）方法开展了基板表面约束阻尼材料的位置和厚度的两级优化。上述文献的研究结果都表明：约束阻尼材料的最优构型与基结构的振动能量存在联系；在优化过程中，都需预先确定材料的几何参数。事实上，几何参数对约束阻尼材料最优构型的影响仍需进一步分析。

为了克服棋盘格式等数值问题，上述优化方法都需引入滤波过程，且优化结果存在中间密度单元。Osher等^[10]最早提出了的水平集法（level set method, LSM），很快被引入到结构拓扑优化领域，其基本思想是利用高一维的水平集函数的零水平集（边界）表示结构拓扑优化的几何边界，通过控制水平集演化速度实现结构形状与拓扑优化。理论上，水平集法可以获得清晰的结构边界信息，也不再需引入滤波过程。Allaire等^[11]将水平集法应用于进行连续结构的拓扑优化，基于形状灵敏度分析中获得水平集的演化速度；Ansari等^[12]以模态损耗因子为目标，将水平集法应用于约束阻尼结构的拓扑优化。上述的传统水平集方法在进行拓扑优化时，需求解复杂的哈密顿–雅可比偏微分方程，在每个迭代步都需采用符号距离函数重新初始化水平集函数，且最优拓扑结构对初始构型依赖性较强。Wang等^[13-14] 提出了参数化水平集方法（parametric level set method, PLSM），将水平集函数解耦为基函数与扩展系数乘积和的形式，避免了求解复杂哈密顿–雅可比偏微分方程。通过经典的最小柔度设计问题说明了参数化水平集法具有收敛速度快、对初始构型依赖弱等优点。Pang等^[15]利用参数化水平集法对自由阻尼结构进行优化设计，算例和实验证明了阻尼材料的优化构型与基板的应变能分布的相似性；Cui等^[16]在利用参数化水平集进行结构优化时，通过增加形状灵敏度约束因子来控制优化步长，提高了结构拓扑优化的效率和稳定性。

基于参数化水平集法的上述优点，本文将参数化水平集法应用于约束阻尼结构的拓扑优化研究，以模态损耗因子最大化为目标函数，约束阻尼材料用量为约束条件，建立约束阻尼结构的拓扑优化模型。基于伴随向量法进行灵敏度分析，采用优化准则法更新设计变量。通过数值算例验证了参数化水平集法的可行性和有效性，着重讨论了约束阻尼材料的初始构型和材料几何参数的变化对最优构型的影响。

1 约束阻尼板有限元模型

采用约束阻尼复合单元建立约束阻尼板的有限元模型^[17]，如图1所示。其中，h表示各层厚度；l和w表示单元节点的坐标符号。每个单元包括4个节点，每个节点有7个自由度，分别是：基板层中性面内位移分量 ${u_{\rm{b}}}$ 和 ${v_{\rm{b}}}$ ；约束层中性面内的两个位移分量 ${u_{\rm{c}}}$ 和 ${v_{\rm{c}}}$ ；约束阻尼板的横向位移 $w$ 以及绕x轴和y轴的转角 ${\varphi _{x}}$ 和 ${\varphi _{y}}$ 。

基于能量法得到基板、约束层和阻尼层的单元质量矩阵和单元刚度矩阵^[17]，并进行组装，那么，约束阻尼板结构的整体质量矩阵和刚度矩阵可表示为

$\left\{ {\begin{array}{*{20}{l}} {{\boldsymbol{M}}={{\boldsymbol{M}}_{\rm{b}}} + \displaystyle\sum\limits_{i=1}^n {\left( {{\boldsymbol{m}}_{{\rm{v}}i} + {\boldsymbol{m}}_{{\rm{c}}i}} \right)} } \\ {{{\boldsymbol{K}}_{\rm{R}}}={{\boldsymbol{K}}_{\rm{b}}} + \displaystyle\sum\limits_{i=1}^n {\left( {{\boldsymbol{k}}_{{\rm{v}}i} + {\boldsymbol{k}}_{{\rm{c}}i} + {\boldsymbol{k}}_{{\rm{\gamma v\_R}}i}} \right)} } \\ {{{\boldsymbol{K}}_{\rm{I}}}=\displaystyle\sum\limits_{i=1}^n {{\boldsymbol{k}}_{{\rm{\gamma v\_I}}i}} } \end{array}} \right.$

(1)

式中： ${\boldsymbol{M}}$ 为约束阻尼板的整体质量矩阵； ${{\boldsymbol{M}}_{\rm{b}}}$ ， ${{\boldsymbol{K}}_{\rm{b}}}$ 为基板层的整体质量矩阵和刚度矩阵； ${\boldsymbol{m}}_{{\rm{v}}i}$ ， ${\boldsymbol{k}}_{{\rm{v}}i}$ ， ${\boldsymbol{m}}_{{\rm{c}}i}$ ， ${\boldsymbol{k}}_{{\rm{c}}i}$ 分别为阻尼层和约束层单元i的质量矩阵和刚度矩阵；n为单元个数； ${\boldsymbol{k}}_{{\rm{\gamma v\_R}}i}$ ， ${\boldsymbol{k}}_{{\rm{\gamma v\_I}}i}$ 分别为阻尼层剪切刚度矩阵的实部和虚部； ${{\boldsymbol{K}}_{\rm{R}}}$ 和 ${{\boldsymbol{K}}_{\rm{I}}}$ 为整体刚度矩阵的实部和虚部。

从而，可以写出约束阻尼板结构的自由振动动力学方程

${\boldsymbol{M}}{\ddot{\boldsymbol{ X}}} + ({{\boldsymbol{K}}_{\rm{R}}} +{\rm{i}}{{\boldsymbol{K}}_{\rm{I}}}){\boldsymbol{X}}=0$

(2)

其中， ${\boldsymbol{X}}$ 为全局节点位移向量。那么，约束阻尼板结构的特征方程为：

$\left( { - \lambda _{r}^{\rm{*}}{\boldsymbol{M}} + ({{\boldsymbol{K}}_{\rm{R}}} + {\rm{i}}{{\boldsymbol{K}}_{\rm{I}}})} \right){\boldsymbol{\varPsi }}_{r}^{\rm{*}}={\boldsymbol{0}}$

(3)

式中， $\lambda _{r}^{\rm{*}}$ 和 ${\boldsymbol{\varPsi }}_{r}^{\rm{*}}$ 分别表示第r阶复特征值和复特征向量。

通常采用模态损耗因子表征约束阻尼结构的抑振能力。定义第r阶的模态损耗因子为

${\eta _r}=\frac{{{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{I}}}{{\boldsymbol{\varPsi }}_{r}}}}{{{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{R}}}{{\boldsymbol{\varPsi }}_{r}}}}$

(4)

式中， ${{\boldsymbol{\varPsi }}_{r}}$ 表示第r阶实模态振型向量。

2 基于参数化水平集法的约束阻尼结构优化设计 2.1 参数化水平集法

水平集法采用水平集函数的零等值线或面表达结构边界，通过追踪水平集函数描述的曲面运动，确定结构边界的演化^[18]。那么，基板表面敷设的约束阻尼材料拓扑形状可以用标准水平集函数 ${\boldsymbol{\varPhi }}({\boldsymbol{x}})$ 进行描述。

$ \left\{ {\begin{array}{*{20}{l}} {{\boldsymbol{\varPhi }}({\boldsymbol{x}}){\rm{ > }}0,\;\;\;{\boldsymbol{x}} \in D{\rm{ }}} \\ {{\boldsymbol{\varPhi }}({\boldsymbol{x}}){\rm{=}}0,\;\;\;{\boldsymbol{x}} \in \varGamma {\rm{ }}} \\ {{\boldsymbol{\varPhi }}({\boldsymbol{x}}){\rm{ < }}0,\;\;\;{\boldsymbol{x}} \in \varOmega /D{\rm{ }}} \end{array}} \right.\;\;\;\; \begin{array}{c}{(}{ ext{材料}}{)}\\ {(}{ ext{边界}}{)}\\ {(}{ ext{孔洞}}{)}\end{array} $

(5)

式中： $\varOmega $ 为设计域； $\varGamma $ 为材料边界； $D$ 为最终的材料区域； ${\boldsymbol{x}}$ 为结构设计域内的节点坐标。水平集函数及材料边界描述如图2。

将优化迭代构型中的隐式水平集函数 ${\boldsymbol{\varPhi }}({\boldsymbol{x}},t)$ 解耦成显式的径向基函数 ${g_{j}}({\boldsymbol{x}})$ 和扩展系数 ${\alpha _{j}}(t)$ 形式方程^[19]，其表达式为

${\boldsymbol{\varPhi }}({\boldsymbol{x}},t)=\sum\limits_{j=1}^m {{g_{j}}\left( {\boldsymbol{x}} \right){\alpha _{j}}(t)} $

(6)

式中：m为设计域内径向基函数个数； $t$ 为伪时间变量。将式（6）写成矩阵形式为

${\boldsymbol{\varPhi }}({\boldsymbol{x}},t)={\boldsymbol{G}}{\rm{(}}{\boldsymbol{x}}{\rm{)}}{\boldsymbol{\alpha }}{\rm{(}}t{\rm{)}}$

(7)

本文采用一种紧支径向基函数^[18] ${g_{j}}\left( {\boldsymbol{x}} \right)$ ，其表达式为

${g_{j}}\left( {\boldsymbol{x}} \right)=\left\{ {\begin{array}{*{20}{l}} {\rm{0}}&{r_{j}}{\rm{(}}{\boldsymbol{x}}{\rm{)}} \geqslant 1 \\ {{\left( {{\rm{1 - }}{r_{j}}{\rm{(}}{\boldsymbol{x}}{\rm{)}}} \right)}^4}\left( {4{r_{j}}{\rm{(}}{\boldsymbol{x}}{\rm{)}} + 1}\right)&{r_{j}}{\rm{(}}{\boldsymbol{x}}{\rm{)}} < 1 \end{array}} \right.$

(8)

式中， ${r_{j}}{\rm{(}}{\boldsymbol{x}}{\rm{)}}$ 为径向基函数的支撑半径，定义为

${r_{j}}{\rm{(}}{\boldsymbol{x}}{\rm{)=}}\frac{{\left\| {{\boldsymbol{x}} - {{\boldsymbol{x}}_{j}}} \right\|}}{R}$

(9)

式中，R为径向基函数的影响半径，取2.5。

2.2 优化模型

在设计域内定义参数化水平集函数后，约束阻尼板结构的质量矩阵和刚度矩阵可以表示为

$\left\{ {\begin{array}{*{20}{l}} {{\boldsymbol{M}}\left( {\boldsymbol{\varPhi }} \right)={{\boldsymbol{M}}_{\rm{b}}} + {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}\left( {{\boldsymbol{m}}_{{\rm{v}}i} + {\boldsymbol{m}}_{{\rm{c}}i}} \right)} \\ {{{\boldsymbol{K}}_{\rm{R}}}\left( {\boldsymbol{\varPhi }} \right)={{\boldsymbol{K}}_{\rm{b}}} + {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}\left( {{\boldsymbol{k}}_{{\rm{v}}i} + {\boldsymbol{k}}_{{\rm{c}}i} + {\boldsymbol{k}}_{{\rm{\gamma v\_R}}i}} \right)} \\ {{{\boldsymbol{K}}_{\rm{I}}}\left( {\boldsymbol{\varPhi }} \right)={\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}{\boldsymbol{k}}_{{\rm{\gamma v\_I}}i}} \end{array}} \right.$

(10)

式中， ${\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}$ 为Heaviside函数，为方便求其偏导数，定义为^[20]

${\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)=}}\left\{ {\begin{array}{*{20}{l}} {a\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {\boldsymbol{\varPhi }} \leqslant - \varDelta } \\ {\dfrac{{3(1 - a)}}{{4\varDelta }}\left( {\dfrac{{\boldsymbol{\varPhi }}}{\varDelta } - \dfrac{{{{\left( {\boldsymbol{\varPhi }} \right)}^3}}}{{3{\varDelta ^3}}}} \right) + \dfrac{{1 + a}}{2}\;\quad - \varDelta < {\boldsymbol{\varPhi }} < \varDelta } \\ {1\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {\boldsymbol{\varPhi }} \geqslant \varDelta } \end{array}} \right.$

(11)

本文中取 $a=0.0001,\varDelta=1$ 。

那么，以参数化水平集函数的扩展系数为设计变量，模态损耗因子最大化为目标函数，建立约束阻尼结构的拓扑优化模型，表述为

$\begin{split} & find:\,\;{\alpha _j},\;\;j=1,2, \cdots, m \\ &\max :\;\,J=\displaystyle\sum\limits_{r=1}^k {{\beta _{r}}{\eta _{r}}} \\ & s.t.\quad {\rm{ }}\left( {{{\boldsymbol{K}}_{\rm{R}}} + {\rm{i}}{{\boldsymbol{K}}_{\rm{I}}}} \right){{\boldsymbol{\varPsi }}_{r}}=\lambda _{r}^{\rm{*}}{\boldsymbol{M}}{{\boldsymbol{\varPsi }}_{r}} \\ & \quad \;\quad {\rm{ }}{\boldsymbol{\varPsi }}_{e}^{\rm{T}}{\boldsymbol{M}}{{\boldsymbol{\varPsi }}_{f}}={\delta _{{ef}}} \\ & \qquad \;{\alpha _{{j\rm{,min}}}} \leqslant {\alpha _{j}} \leqslant {\alpha _{{j\rm{,max}}}} \\ &\qquad\; V{\rm{=}}{\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}\sum\limits_{i=1}^n {{v_{i}}} \leqslant \bar V \end{split} $

(12)

式中： $V$ 为基板表面敷设的约束阻尼材料体积； $\bar V$ 为体积约束； ${v_{i}}$ 为单个单元i内约束阻尼材料的体积；e，f均为阶数； $ \delta $ 是对 ${H}{\rm{(}}{\varPhi } {\rm{)}}$ 进行光滑化处理的区域界限； ${\alpha }_{j,\rm{min}} ext{、}{\alpha }_{j,\rm{max}}$ 为设计变量的上下限，取 ${\alpha _{{j,\rm{min}}}}=2\min \alpha ,{\alpha _{{j,\rm{max}}}}=2\max \alpha$ ； ${\ \beta _{r}}$ 是权重系数。

2.3 灵敏度分析

本文采用伴随向量法进行灵敏度分析。为避免计算目标函数中振型对设计变量的偏导数，以及方便采用优化准则法更新设计变量，通过引入伴随向量 ${{\boldsymbol{\ \mu }}_{\rm{1}}}$ 和 ${\ \mu _2}$ 将目标函数改写为

$\begin{split}{J_{r}}=& \frac{{{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{R}}}{{\boldsymbol{\varPsi }}_{r}}}}{{{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{I}}}{{\boldsymbol{\varPsi }}_{r}}}}{\rm{ + }}{\rm{R_{\rm{e}}}}\left( {{\boldsymbol{\mu }}_{\rm{1}}^{\rm{H}}\left( {{\boldsymbol{K}}_{\rm{R}}+{\rm i}{\boldsymbol{K}}_{\rm{I}} - \lambda _{r}^{\rm{*}}{\boldsymbol{M}}} \right){{\boldsymbol{\varPsi }}_{r}}} \right) + \\ &{\mu _{\rm{2}}}\left( {{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{\boldsymbol{M}}{{\boldsymbol{\varPsi }}_{r}} - 1} \right)\end{split}$

(13)

式中， ${\rm{R_{\rm{e}}}}$ 表示取实部。

那么，目标函数对设计变量 ${\alpha _{j}}$ 的灵敏度可表示为

$\begin{split}\dfrac{{\partial {J_{r}}}}{{\partial {\alpha _{j}}}}=&\dfrac{{\left\{ \begin{array}{l} \left( {{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{I}}}{{\boldsymbol{\varPsi }}_{r}}} \right)\left( {{\boldsymbol{\varPsi }}_{r}^{\rm{T}}\dfrac{{\partial {{\boldsymbol{K}}_{\rm{R}}}}}{{\partial {\alpha _{j}}}}{{\boldsymbol{\varPsi }}_{r}}} \right) \\ - \left( {{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{R}}}{{\boldsymbol{\varPsi }}_{r}}} \right)\left( {{\boldsymbol{\varPsi }}_{r}^{\rm{T}}\dfrac{{\partial {{\boldsymbol{K}}_{\rm{I}}}}}{{\partial {\alpha _{j}}}}{{\boldsymbol{\varPsi }}_{r}}} \right) \end{array} \right\}}}{{{{\left( {\phi _{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{I}}}{\phi _{r}}} \right)}^2}}}{\rm{ + }}\\ &{\rm{R_{\rm{e}}}}\left( {{\boldsymbol{\mu }}_{\rm{1}}^{\rm{H}}\left( {\dfrac{{\partial {\boldsymbol{K}}}}{{\partial {\alpha _j}}} - \lambda _{r}^{\rm{2}}\dfrac{{\partial {\boldsymbol{M}}}}{{\partial {\alpha _{j}}}}} \right){{\boldsymbol{\varPsi }}_{r}}} \right) + {\mu _2}\left( {{\boldsymbol{\varPsi }}_{r}^{\rm{T}}\dfrac{{\partial {\boldsymbol{M}}}}{{\partial {\alpha _{i}}}}{{\boldsymbol{\varPsi }}_{r}}} \right)\end{split}$

(14)

式中，伴随向量 ${{\boldsymbol{\mu }}_{\rm{1}}}$ 、 ${\mu _2}$ 是下列方程的解。

$\begin{split} &\left[ {\begin{array}{*{20}{c}} {\left( {{\boldsymbol{K}} - \lambda _{r}^{\rm{2}}{\boldsymbol{M}}} \right)}&{2{\boldsymbol{M}}{{\boldsymbol{\varPsi }}_{r}}} \\ {2{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{\boldsymbol{M}}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\mu }}_1}} \\ {{\mu _2}} \end{array}} \right]=\\ &\qquad\left[{ \begin{array}{*{20}{c}} \dfrac{ - \left( {{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{I}}}{{\boldsymbol{\varPsi }}_{r}}} \right)\left( {2{{\boldsymbol{K}}_{\rm{R}}}{{\boldsymbol{\varPsi }}_{r}}} \right) + \left( {{\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{I}}}{{\boldsymbol{\varPsi }}_{r}}} \right)\left( {2{{\boldsymbol{K}}_{\rm{R}}}{{\boldsymbol{\varPsi }}_{r}}} \right)}{\left( {\boldsymbol{\varPsi }}_{r}^{\rm{T}}{{\boldsymbol{K}}_{\rm{I}}}{{\boldsymbol{\varPsi }}_{r}} \right)^2 }\\ 0 \end{array}} \right] \end{split}$

(15)

质量矩阵和刚度矩阵对设计变量的灵敏度为

$\left\{ {\begin{array}{*{20}{l}} {\dfrac{{\partial {\boldsymbol{K}}\left( {\boldsymbol{\varPhi }} \right)}}{{\partial {\alpha _{j}}}}{\rm{=}}\dfrac{{\partial {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}}}{{\partial {\boldsymbol{\varPhi }}}}\dfrac{{\partial {\boldsymbol{\varPhi }}}}{{\partial {\alpha _{j}}}}\left( {{\boldsymbol{k}}_{{\rm{c}}i} + {\boldsymbol{k}}_{{\rm{v}}i} + {\boldsymbol{k}}_{{\rm{\gamma v\_R}}i}{\rm{ + }}{\rm{i}}{\boldsymbol{k}}_{{\rm{\gamma v\_I}}i}} \right){\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}} \\ {\dfrac{{\partial {\boldsymbol{M}}\left( {\boldsymbol{\varPhi }} \right)}}{{\partial {\alpha _{j}}}}{\rm{=}}\dfrac{{\partial {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}}}{{\partial {\boldsymbol{\varPhi }}}}\dfrac{{\partial {\boldsymbol{\varPhi }}}}{{\partial {\alpha _{j}}}}\left( {{\boldsymbol{m}}_{{\rm{v}}i} + {\boldsymbol{m}}_{{\rm{c}}i}} \right){\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}} \\ {\dfrac{{\partial {{\boldsymbol{K}}_{\rm{R}}}\left( {\boldsymbol{\varPhi }} \right)}}{{\partial {\alpha _{j}}}}{\rm{=}}\dfrac{{\partial {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}}}{{\partial {\boldsymbol{\varPhi }}}}\dfrac{{\partial {\boldsymbol{\varPhi }}}}{{\partial {\alpha _{j}}}}\left( {{\boldsymbol{k}}_{{\rm{c}}i} + {\boldsymbol{k}}_{{\rm{v}}i} + {\boldsymbol{k}}_{{\rm{\gamma v\_R}}i}} \right){\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}} \\ {\dfrac{{\partial {{\boldsymbol{K}}_{\rm{I}}}\left( {\boldsymbol{\varPhi }} \right)}}{{\partial {\alpha _{j}}}}{\rm{=}}\dfrac{{\partial {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}}}{{\partial {\boldsymbol{\varPhi }}}}\dfrac{{\partial {\boldsymbol{\varPhi }}}}{{\partial {\alpha _{j}}}}{\boldsymbol{k}}_{{\rm{\gamma v\_I}}i}{\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}} \end{array}} \right.$

(16)

Heaviside函数的偏导数由式（11）可得

$\frac{{\partial {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}}}{{\partial {\boldsymbol{\varPhi }}}}{\rm{=}}\left\{ {\begin{array}{*{20}{c}} {\dfrac{{{\rm{3}}\left( {{\rm{1 - }}a} \right)}}{{{\rm{4}}\varDelta }}\left( {{\rm{1 - }}{{\left( {\dfrac{{\boldsymbol{\varPhi }}}{\varDelta }} \right)}^{\rm{2}}}} \right)\quad \quad \left| {\boldsymbol{\varPhi }} \right| \leqslant \varDelta } \\ {b\quad \quad \quad \quad \quad \quad \quad\quad \quad \left| {\boldsymbol{\varPhi }} \right| > \varDelta } \end{array}} \right.$

(17)

式中： $a=0.000\, 1$ ； $b=0.000\, 1$ ； $ \varDelta=1 $ 。水平集函数对设计变量的导数由式（7）可得

$\frac{{\partial {\boldsymbol{\varPhi }}}}{{\partial {\alpha _{j}}}}{\rm{=}}{\left[ {\begin{array}{*{20}{c}} {g_1\left( {{{\boldsymbol{x}}}} \right)}& \cdots &{g_m\left( {{{\boldsymbol{x}}}} \right)} \end{array}} \right]^{\rm{T}}}$

(18)

此外，体积约束对设计变量的灵敏度为

$\frac{{\partial V}}{{\partial {{\rm{\alpha }}_{j}}}}=\frac{{\partial {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}}}{{\partial {{\rm{\alpha }}_{j}}}}\sum\limits_{i=1}^n {{v_{i}}}=\frac{{\partial {\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}}}{{\partial {\boldsymbol{\varPhi }}}}\frac{{\partial {\boldsymbol{\varPhi }}}}{{\partial {{\rm{\alpha }}_{j}}}}\sum\limits_{i=1}^n {{v_i}} $

(19)

2.4 设计变量更新准则

优化准则（optimization criterion，OC）法在求解带约束条件的优化问题时，先构造拉格朗日函数，再利用库恩-塔克条件（KKT）构造极值条件函数，进而确定设计变量更新准则，构造约束条件和目标函数的拉格朗日函数

$\begin{split}L=& {J_{r}}\left( {{\alpha _{j}}} \right) + \lambda _{j}^{\rm{1}}\left({\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)}}\sum\limits_{i=1}^n {{v_{i}}} - \bar V\right) + \\ &\sum\limits_{j=1}^m {\lambda _{j}^{\rm{2}}({\alpha _{j,{\rm{min}}}} - {\alpha _{j}})} + \sum\limits_{j=1}^m {\lambda _{j}^{\rm{3}}({\alpha _{j}} - {\alpha _{j,{\rm{max}}}})} \end{split}$

(20)

式中, $\lambda _{j}^{\rm{1}}$ 、 $\lambda _{j}^{\rm{2}}$ 、 $\lambda _{j}^{\rm{3}}$ 为拉格朗日乘子。由KKT条件，确定的设计变量更新准则为

$\alpha _{j}^{{(q + 1)}}{\rm{=}}\left\{ {\begin{array}{*{20}{l}} {{\rm{max(0,}}\;{\alpha _{j}} - p{\rm{), }}\;\;\;\;\;{\alpha _{j}} < {\alpha _{{\rm{min}}}}} \\ {{{\left( {D_{j}^{{(q)}}} \right)}^{\rm{\xi }}}\alpha _{j}^{{(k)}}{\rm{, }}\;\;\;\;\;\;\;{\alpha _{{\rm{min}}}} \leqslant {\alpha _{j}} \leqslant {\alpha _{{\rm{max}}}}} \\ {{\rm{min(1,}}\;{\alpha _{j}} + p{\rm{), }}\;\;\;\;{\alpha _{i}} > {\alpha _{{\rm{max}}}}} \end{array}} \right.$

(21)

式中： $D_{j}^{{(q)}}{\rm{=}}{{\dfrac{{\partial {J_{r}}({\alpha _{j}})}}{{\partial {\alpha _{j}}}}} / {{\lambda _{j}}}}V$ ；q为迭代次数。为了保证设计变量迭代过程稳定，引入两个参数，分别是阻尼因子 $\xi (0 < \xi < 1)$ 和移动极限 $p$ ， $p$ 用来控制迭代的步长， $\xi $ 取0.3，前50次迭代 $p$ 值取0.015，超过50次取0.001 5。

3 优化流程

基于参数化水平集法的约束阻尼板壳结构的拓扑优化流程如图3所示，具体描述如下：

a. 初始化水平集函数 $\varPhi$ 、径向基函数 ${\boldsymbol G}\left( {\boldsymbol{x}} \right)$ ，并依据式（7）初始化设计变量（扩展系数） ${{{\alpha }}_{j}}$ ；

b. 依据初始化的水平集函数确定约束阻尼材料的初始构型，并建立对应的约束阻尼板结构的有限元模型；求解约束阻尼结构的特征方程，得到初始结构的模态损耗因子；

c.采用伴随向量法分析目标函数对设计变量的灵敏度；

d. 依据优化准则法更新设计变量，即参数化水平函数的扩展系数，得到更新的水平集函数，更新约束阻尼板结构的有限元模型和特征方程；

e.求解特征方程得到更新的目标函数以及体积约束条件；

f. 判断收敛条件是否满足，若不满足，重复步骤c～e，直至满足迭代终止条件。本文中的收敛条件为：体积分数与目标体积之差、目标函数值连续20代的变化值同时小于0.001。

4 算例分析和讨论

基于上述的建模方法以及参数化水平集法，编制了约束阻尼板结构的有限元模型及拓扑优化程序。约束阻尼板的参数如表1所示，基板尺寸为0.4 m $ imes $ 0.3 m。粘弹性阻尼材料的剪切模量采用复常数模型，损耗因子取0.5，边界条件为四边固支，体积约束为全覆盖约束阻尼材料的50%，有限元单元数为 $40 imes 30$ 。分别以一阶和二阶模态损耗因子及其加权最大化为优化目标，采用参数化水平集法对基板表面的约束阻尼材料进行拓扑优化。

4.1 优化结果与分析

已有文献研究表明，以模态损耗因子为目标函数对约束阻尼结构进行拓扑优化时，约束阻尼材料优化构型与基板的模态应变能分布有一定的相似性^[15]。本文中基板的前二阶模态应变能及归一化加权模态应变能分布如图4。

此外，基于参数化水平集法的结构拓扑优化对初始结构依赖性较弱^[14]，本文以一阶模态损耗因子最大化进行拓扑优化时，构建了约束阻尼材料覆盖率约为50%（体积约束）的两种初始构型：a. 基板表面均匀覆盖若干约束阻尼材料，记为初始构型A，如图5（a）所示；b. 约束阻尼材料依据基板单元应变能分布构建初始构型，记为初始构型B，如图5（b）所示。

采用上述两种初始构型，优化得到的水平集函数相同，如图6（a）。进一步，将图6（a）的水平集函数规整后获得约束阻尼材料在基板表面的最优布局，如图6（b）所示（蓝色表示约束阻尼材料，白色表示基板结构）。本文中水平集函数规整准则为：

${\boldsymbol{H}}{\rm{(}}{\boldsymbol{\varPhi }}{\rm{)=}}\left\{ {\begin{array}{*{20}{c}} {0,\;\;\;{\boldsymbol{\varPhi }} \leqslant {{d}}} \\ {1,\;\;\;{\boldsymbol{\varPhi }} > {{d}}} \end{array}} \right.$

(22)

式中, ${d}$ 表示对水平集函数排序后第 $\bar V imes {n_{{\rm{nel}}}}$ 单元的水平集函数值， ${n_{{\rm{nel}}}}$ 是总单元数。显然，基于参数化水平集法的约束阻尼拓扑优化对初始构型的选取依赖性很弱。图7所示为两种初始构型的优化历程，初始构型A优化用时233 s，初始构型B优化用时188 s，可以看出，依据模态应变能构建初始构型，可以提高拓扑优化效率。因此，下文中约束阻尼材料的初始布局都基于基板结构的模态应变能分布构建。

为表征最优构型与依据基板应变能构建的初始构型相关性，本文定义：约束阻尼材料最优构型和初始构型的重合面积与初始构型的面积之比为相似度。图6的最优构型与初始构型相似度为0.7933，表明在表1参数下，约束阻尼材料的最优布局与初始结构布局B的相似度较高。但对比图5（b）和图6（b）可知，在应变能比较大的板结构中心区域，约束阻尼材料敷设面积明显减小，优化后的一阶模态损耗因子为0.0259，与初始结构布局B对应的模态损耗因子相比较，提高了14.1%。

以二阶模态损耗因子最大化为目标函数时，基于基板第二阶模态应变能分布构建初始布局，如图8（a），相应的约束阻尼材料最优布局如图8（b），优化迭代过程如图8（c）。图8中优化构型与初始构型相似度为0.84。对比图8（a）和（b），同样可以看出，在应变能比较大的板结构中心区域，约束阻尼材料敷设面积明显减小，优化后的二阶模态损耗因子为0.0256，比初始结构布局对应的模态损耗因子提高了8.1%。

进一步，以一阶和二阶加权模态损耗因子最大化为目标函数，基于基板的一阶和二阶加权模态应变能分布构建初始布局，如图9（a）。利用参数化水平集法进行加权优化，约束阻尼材料最优布局如图9（b），优化迭代过程如图9（c）。权重系数取 ${\ \beta _1}{\rm{=}}{\beta _2}{\rm{=}}0.5$ 。

图9中优化构型与初始构型相似度为0.8，优化后的加权模态损耗因子为0.024，比初始结构布局对应的模态损耗因子提高了14%。从优化历程也可看出，以基板模态应变能构建初始构型时进行的拓扑优化效率较高。表2列出了上述算例中的模态损耗因子优化结果，可以看出，与初始构型相比较，单一目标下的模态损耗因子都有所提高；考虑加权目标函数时，尽管第二阶的模态损耗因子较初始结构有所降低，但仍能很好地兼顾前两阶模态的减振性能。综上所述，基于参数化水平集法的约束阻尼结构拓扑优化方法有效可行。

4.2 优化结果对厚度参数的依赖性分析

进一步讨论在附加质量保持不变的条件下，约束层材料和阻尼材料的厚度变化对约束阻尼材料模态损耗因子和最优构型的影响规律。

保持约束阻尼材料的附加质量不变，即体积分数（质量分数）为0.5，改变约束层和阻尼层的厚度，以一阶和二阶及其加权模态损耗因子为目标函数进行优化，分析拓扑优化构型以及模态损耗因子变化。优化结果对应的模态损耗因子和相似度变化如图10。

从图10可以看出，在保持约束阻尼层附加质量相同的条件下，阻尼层厚度与约束层厚度的比值逐渐减小时，优化构型的相似度也减小，即最优构型与模态应变能分布的布局差异变大。但是存在一个最优的约束层和阻尼层厚度比，使得模态损耗因子最大，此时优化构型的相似度都大于0.75，且约束阻尼材料在基板中心区域覆盖面积减小，四边覆盖面积增加。因此，在上述附加质量的设计约束下，当约束层和阻尼层厚度接近相同时，可依据上述约束阻尼材料分布特征构建初步构型。

保持约束阻尼层附加质量不变，同时设约束阻尼材料的总厚度不变，改变约束层和阻尼层的相对厚度。这种情况下，虽然附加质量不变，但附加的体积分数不同，如表3。以一阶和二阶及其加权模态损耗因子为目标函数进行优化，不同约束层和阻尼层厚度对应的模态损耗因子和相似度变化如图11。

从图11可以看出，在保持约束阻尼层的厚度和附加质量不变的条件下，随着阻尼层与约束层厚度之比的逐渐减小，模态损耗因子先增大后减小，存在一个最适厚度比使得模态损耗因子最大，在本文的参数下，约束层与阻尼层最适厚度比约为7∶1。当选取最适厚度时，优化构型的相似度大于0.75。这为在附加质量和材料总厚度的设计约束下，提供了约束阻尼材料设计的参考思路。

5 结　论

结合约束阻尼结构的有限元模型，建立了基于参数化水平集法的约束阻尼结构动力学拓扑优化模型，以前两阶及其加权模态损耗因子最大化为目标函数进行了拓扑优化。结果表明：基于参数化水平集法对约束层阻尼结构进行拓扑优化是可行有效的，且以基板结构应变能分布构建初始构型可提高优化效率。

在保持约束阻尼材料附加质量相同的条件下，阻尼层厚度与约束层厚度的比值逐渐减小时，优化构型的相似度也减小；但是存在一个最优的约束层和阻尼层厚度比，使得模态损耗因子最大，优化构型的相似度大于0.75。在保持约束阻尼层的厚度和附加质量不变的条件下，随着阻尼层与约束层厚度之比逐渐减小时，模态损耗因子先增大后减小，存在一个最适厚度比使得模态损耗因子最大，在本文的参数下，约束层与阻尼层最适厚度比约为7:1。当选取最适厚度时，优化构型的相似度大于0.75。这为在附加质量和材料总厚度的设计约束下，提供了约束阻尼材料设计的参考思路。