水下武器和作战平台的周期性机械振动会使结构与声介质产生耦合作用而向外辐射噪声，影响武器系统的隐身性和自身的声呐探测距离．约束阻尼层在减振降噪效果和有效作用频带等方面有优势，在振动和噪声的控制方面得到广泛应用，但增加了结构质量．如何在尽量减少阻尼材料用量的前提下降低振动噪声成为了关键，因此结合连续体拓扑优化法合理的设置阻尼层分布以达到最佳的降噪效果具有重要意义．拓扑优化法用于降低结构噪声，主要分为基于声振耦合有限元法的内声场优化和基于有限元法/边界元法的辐射噪声优化．在优化内声场的研究中，文献[1-4]以材料体积为约束条件，以目标区域声压最小化为目标函数进行优化设计．由于有限元法的局限性，仅能处理有限空间的优化设计问题．对于辐射声场优化问题一般结合有限元法和边界元法，以空气为声介质，以阻尼材料的相对密度为设计变量，使辐射声功率或空间某点声压最小化．文献[5]根据经典薄板理论，建立约束阻尼板有限元模型，利用瑞利积分推导辐射声功率及灵敏度表达式．文献[6]同样以瑞利积分推导辐射声功率，对比模态损耗因子最大化和声压最小化为目标的附加自由阻尼板优化结果，认为以声压为目标的优化效果更好．文献[7]以结构表面辐射声功率最小为设计目标，研究复合材料安静结构设计方法问题．文献[8]在对吸声材料层优化声压的研究中，发现吸声层的能量损失效应并提出声动柔度的概念．文献[9]研究了非比例阻尼作用下结构声学性能的灵敏度分析及壳结构阻尼层的优化设计，并提出了一种用于声压灵敏度分析的伴随变量法．空气中结构声辐射可以忽略声介质对板的耦合作用，从而声功率和声压灵敏度容易通过解耦求取．但是工程中还存在潜艇等水下机械结构，水对结构的耦合作用是不可忽略的，采用传统的弱耦合模型难以满足水下结构的声学设计需求．这里建立了强耦合含弹性边界的声振耦合有限元/边界元模型，提出了一种强耦合声功率和声压灵敏度求解方法．此外，根据质量和刚度矩阵与薄板厚度的关系，建立以薄板结构上阻尼材料厚度为设计变量的人工阻尼模型．通过不同设计条件下的声辐射优化设计结果验证了灵敏度方法和优化方法的适用性．1 结构振动/声辐射方程无限大空间中的边界积分方程Hpb-Gvn=0，(1)式中：H和G为采用边界积分方程得到的频率相关的矩阵；pb和vn为边界上的声压和法向速度矩阵．在外载荷角频率为ω的谐激励作用下，强耦合问题须要引入声介质对结构表面的反作用项，薄板阻尼结构的振动方程满足(Ks-iωC-ω2Ms)Xs=Fs+Lpf，(2)式中：Ks为结构刚度矩阵；C为结构的阻尼矩阵；Ms为结构质量矩阵；Xs为结构位移响应向量；Fs为结构所受外载荷；L为耦合面上的力耦合矩阵；pf为场点声压向量；i为虚部符号．安装有弹性减震器的板结构的总刚度矩阵由薄板结构的连续单元刚度矩阵和弹簧单元的刚度矩阵组成．弹簧质量较小，若忽略弹簧质量对结构总质量矩阵的影响，则结构系统的质量矩阵和刚度矩阵可表示为Ms=∑i=1nmshi;Ks=∑i=1nkshi+∑j=1mkspj,式中：mshi和kshi分别为第i个板单元的单元质量矩阵和单元刚度矩阵；n为板单元总数；kspj为第j个弹簧单元的刚度矩阵；m为弹簧总数．声介质的法向速度即结构位移在外法向的分量对时间t的导数，则耦合面上声介质法向速度与结构位移向量的关系为[10]vn=-iωD-1LTXs，(3)式中D=∫SNfTNfdS为耦合边界矩阵，Nf为单元形函数．阻尼板与声介质的声振耦合方程为Ks-iωC-ω2Ms-L       iωGD-1LTHXspb=Fs0．(4)通过式(4)可求出耦合面上的声压和薄板位移响应，得到场内任一点的声压pf和薄板的辐射声功率W，具体为：pf=Gqvn-Hqpb；(5)W=-Re(pbTDvnH)/2，(6)式中Gq和Hq为场点中的边界元系数矩阵．为方便求单元灵敏度，将式(1)代入式(6)并由复共轭矩阵乘法规则得到采用边界振速表示的辐射声功率为W=-(vnHBvn)/4，(7)式中B=(H-1G)HD+((H-1G)HD)H．将式(3)代入式(7)，得到辐射声功率关于Xs的表达形式为W=-(XsHQXs)/4，(8)式中Q=(iωD-1LT)HB(iωD-1LT)．2 单元灵敏度分析以阻尼层单元的相对厚度xi为设计变量时须确定设计变量与单元质量、刚度和阻尼矩阵的关系．Kirchhoff板单元的质量与厚度ti呈正比，刚度与厚度ti的三次方呈正比[11]，即mid/m0d=ti/t0=xi;kid/k0d=(ti/t0)3=xi3,式中：t0为最大厚度；mid和kid分别为第i个阻尼单元的单元质量矩阵和单元刚度矩阵；m0d和k0d分别为初始阻尼单元的单元质量矩阵和单元刚度矩阵．含阻尼层的薄板单元的单元质量矩阵me和刚度矩阵ke可视为阻尼单元和基体层单元的叠加，有：me=m0b+xim0d;ke=k0b+xi3k0d. (9)基体层的阻尼与阻尼层相比很小，在此忽略基体层的阻尼，并假设阻尼层的阻尼符合比例阻尼特性．根据式(9)质量和刚度与相对厚度的关系可得Cie=xiα0m0d+xi3β0k0d，(10)式中：Cie为第i个单元的单元阻尼矩阵；α0和β0为瑞利阻尼系数．2.1　辐射声功率灵敏度分析声功率灵敏度是声功率对xi求微分．式(8)两边对xi求导可得∂W/∂xi=-(QXs∂XsH/∂xi+XsHQ∂Xs/∂xi)/4=-Re(XsHQ∂Xs/∂xi)/2. (11)式(1)、式(2)和式(3)分别对xi求导可得：H∂pb/∂xi-G∂vn/∂xi=0;Xs∂Kd/∂xi+Kd∂Xs/∂xi=L∂pb/∂xi;∂vn/∂xi=-iωD-1LT∂Xs/∂xi, (12)式中Kd=Ks-iωC-ω2Ms．求解方程组(12)，经解耦得到∂Xs/∂xi=-(Kd+iωLH-1GD-1LT)-1Xs∂Kd/∂xi ．(13)结合式(9)和(10)可求出Kd对xi的微分为∂Kd/∂xi=3xi2k0d-ω2m0d+iω(α0m0d+3xi2β0k0d). (14)将式(13)代入式(10)和(11)，并结合式(14)得到辐射声功率的单元灵敏度公式为∂W∂xi=12ReXsHQ(Kd+iωLH-1GD-1LT)-1⋅[3xi2k0d-ω2m0d+iω(α0m0d+3xi2β0k0d)]Xs.2.2　声压灵敏度分析结合式(3)和式(5)得到场内一点的声压表达式为pf=iω(-Gq+HqH-1G)D-1LTXs，(15)式(15)中仅Xs随xi发生改变，对xi求微分得到∂pf/∂xi=iω(-Gq+HqH-1G)D-1LT∂Xs/∂xi．(16)将式(13)代入式(16)并结合式(14)得到pf的单元灵敏度表达式为∂pf/∂xi=-iω(-Gq+HqH-1G)D-1LT⋅(Kd+iωLH-1GD-1LT)-1[3xi2k0d-ω2m0d+iω(α0m0d+3xi2β0k0d)]Xs.3 数值算例算例均采用自编程的方式计算出单元质量和单元刚度矩阵，并组装成系统矩阵及实现后续的计算和优化工作．算例1以空气为声介质的强/弱耦合优化对比以空气为声介质的辐射优化问题可以近似地忽略空气对薄板结构的反作用，只考虑弱耦合作用．文献[9]和[12]均采用这种简化模型对阻尼板或多材料的复合板结构的声辐射特性进行了优化．首先以文献[9]提出的弱耦合条件下的伴随变量法为对比算例，验证本研究耦合模型和灵敏度的正确性．图1为阻尼薄板结构耦合系统的示意图，结构部分由长宽为2.0 m、厚度均为0.02 m的基体层和阻尼层组成．板的四个顶角以弹簧接地，弹簧刚度为80 kN/m．声介质为空气，声速为343 m/s，密度为1.21 kg/m3．基体层弹性模量为6.9×1010 Pa，密度为2 700 kg/m3，泊松比为0.3，阻尼层弹性模量为2.35×109 Pa，密度为1 180 kg/m3，泊松比为0.38．本算例以板中点上方3 m处的声压值最小化为目标函数，设置目标体积为全结构域的50%不变．计算域用0.05 m的四边形板单元离散，过滤半径rmin=0.1 m，最小相对高度xmin=1×10-4，单元状态变化量为Δxi=0.03．10.13245/j.hust.210921.F001图1含弹性边界的阻尼薄板结构耦合系统选取阻尼参数α0=2.3，β0=0.02，阻尼系数与频率的关系如图2所示．外载荷垂直于板中点，幅值为1 000 N，频率选取声压响应峰值处的激励频率ω=31 Hz，不同耦合模型优化后的阻尼分布如图3所示．图3(a)采用弱耦合伴随变量法，图3(b)为强耦合条件下本研究推导的灵敏度的优化结果．可以得知两种方法优化后的阻尼分布基本一致．图4为两种方法的迭代过程．强耦合方法的目标函数值略低于弱耦合的响应值，这与理论相符合．因为强耦合考虑了介质对板的反作用，一定程度上抑制了板的振动，所以声压略有降低．数值算法在20个迭代步内达到收敛，说明优化算法的数值稳定性较好．为进一步说明两种耦合方式的降噪效果，表1给出了优化前后参考点声压级．两种耦合模式的降噪幅值均在20 dB左右，效果显著．通过以上对比说明本研究的强耦合模型和声压灵敏度推导正确．10.13245/j.hust.210921.F002图2阻尼系数曲线10.13245/j.hust.210921.F003图3不同耦合模型优化后的阻尼分布10.13245/j.hust.210921.F004图4强/弱耦合目标值迭代过程10.13245/j.hust.210921.T001表1强/弱耦合声压优化效果对比耦合程度初始优化降噪幅值强耦合108.889.619.2弱耦合111.889.822.0dB选取板上方2 m×2 m区域的自由空间内的声压分布情况来对比不同设计的降噪效果，如图5所示，三幅图在同一色度范围内显示以方便对比声压级的大小．从图中虚线标注的高声压外轮廓高度可以看出，弱耦合和强耦合优化后的外场声压级明显减小．两种方法的声压级分布相近，但弱耦合的高声压级区域较强耦合的略大，说明以空气为声介质的阻尼优化问题可以采用弱耦合模型，同时不失优化的准确性．10.13245/j.hust.210921.F005图5不同设计的的外场声压级分布算例2以声压或声功率为目标函数为了探究相同边界条件下以声功率最小化和以目标点声压最小化为目标的函数对阻尼分布的影响，以基体层和阻尼层厚度为0.02 m，尺寸为2 m×2 m的含阻尼覆盖层的简支板为研究对象．声介质为水，密度为1 000 kg/m3，声速为1 450 m/s．基体板材料的弹性模量为2.1×1011 Pa，密度为7 800 kg/m3，阻尼材料的弹性模量为2.36×1010 Pa，密度为3 100 kg/m3，泊松比均为0.3．选取阻尼系数α0=1.2，β0=0.01，外载荷垂直于板中点，幅值为1 000 N，激励频率ω=79 Hz，设置目标体积为全域的50%不变．过滤半径rmin=0.08 m，最小相对高度xmin=1×10-4，单元状态变化量为Δxi=0.03．设置7个以不同目标点声压值最小化的算例，另一个算例以声功率最小化为目标，目标点位置如图6所示．10.13245/j.hust.210921.F006图6声压优化目标点位置示意图图7为8个算例优化后的阻尼分布，通过对比板中点正上方不同高度的声压优化结果(a)、(b)、(c)与声功率的优化分布(h)可以发现4个算例的阻尼材料分布基本一致．这说明对于规则板，当目标点位于板的几何中心垂直线上时，对优化结果没有影响．观察算例(d)、(e)和(f)可以发现：阻尼层的分布受目标点与中心垂线的偏移量影响，阻尼层分布呈现微弱的左右不对称，且偏向目标点一侧的阻尼分布较多，阻尼分布的非对称性随着偏移量的增大变得更明显．如果偏移量过大(图7(g))，那么阻尼分布的非对称性反而降低．在实际应用中，目标点往往相对机械结构的位置是非唯一的．通过以上对比可以发现：声压优化和声功率优化结果相近，且目标点偏移对阻尼分布并没有太大影响，因此采用以声功率为目标的优化更有实际意义．10.13245/j.hust.210921.F007图7声功率优化和不同目标点声压优化的阻尼分布图8为板上方2 m×2 m区域不同设计的外声场声压级分布云图．由两种不同目标函数优化结果的外场声压级分布可知，声压和声功率最小化的优化设计的外场声压级分布基本相同，与图7的优化结果相符．同时，优化设计实现外声场的高声压的分布范围和图中显示的声压级最高值降低．图9为三种优化设计的声功率级扫频曲线．以声功率为目标的优化扫频与以声压为目标的优化扫频结果完全一致．对比优化前后的声功率级扫频曲线可见，通过合理分布阻尼材料的位置和厚度，优化后的结构很好地抹平了响应的共振峰，实现了共振频率附近响应的降低．10.13245/j.hust.210921.F008图8不同设计的外场声压级分布10.13245/j.hust.210921.F009图9不同设计的声功率级扫频曲线算例3弹簧刚度、激励频率和阻尼系数等因素对声辐射的影响首先修改边界上位移约束弹簧的刚度考察边界支撑状况对辐射声功率的影响，算例模型的尺寸和材料参数和3.2节相同．板的四条边上5个位移约束弹簧均匀分布，弹簧刚度分别取1×105，1×106，1×107，1×108和1×109 N/m，另外设置四边简支的算例为对比算例．垂直于板中点施加外载荷，幅值为1 000 N．图10为不同边界条件下板的声功率级扫频曲线．可以发现：随着弹簧刚度的增加，曲线整体起伏的形态相近，但是响应曲线的峰值频率提高，曲线逐渐逼近简支条件的响应曲线．这说明位移约束弹簧的刚度明显改变板的辐射频率，且辐射频率随刚度增大而增大，直至最后逼近简支边界．10.13245/j.hust.210921.F010图10不同边界下板的声功率级扫频曲线其次研究不同激励频率对优化结果的影响，分别考虑3种不同的激励频率，即ω分别取8，56和145 Hz，以声功率最小化为目标．为提高优化效果，方板采用60×60的矩形Kirchhoff板单元离散，初始阻尼体积分数为0.5，均匀分布．选取弹簧边界的刚度为1×107 N/m，矩形板每边5个位移弹簧单元均匀分布，垂直于板中点施加外载荷．过滤半径取2.5倍单元尺寸，最小相对高度xmin=1×10-4，单元状态变化量为Δxi=0.05．不同激励频率下板结构的振动形态和优化阻尼分布如图11所示，板的振动取垂直于板方向的振动幅值．优化结果显示阻尼层分布与激励频率直接相关．随着激励频率的提高，阻尼分布形态的复杂度也越高，并且阻尼多是覆盖在振幅较高的区域．这是因为激励力激起的振型随着频率的提高也变复杂，阻尼分布须要适应振动高振幅区域以抑制振动，实现降低声辐射．10.13245/j.hust.210921.F011图11不同激励频率下阻尼材料分布进一步探究阻尼系数对降噪效果和优化结果的影响，为控制变量，保持α0不变，通过改变对结构振动影响较大的β0考察阻尼系数对优化结果的影响．设定α0=1.2，β0分别为0.005，0.010，0.020和0.030，激励频率ω=145 Hz，等体积约束，以声功率最小化为目标，优化前后声功率级变化情况如表2所示．随着β0的增大，优化后降噪幅值呈现先增大后减小的趋势，说明并不一定β0越大降噪效果越好，而是存在一个合理的取值范围．图12为不同阻尼系数下阻尼材料的优化分布，可以发现系数β0越大，优化后的阻尼材料越集中．10.13245/j.hust.210921.T002表2不同阻尼系数的优化效果对比阻尼系数声功率级/dB降噪幅值/dB初始优化α0=1.2，β0=0.010123.0110.612.4α0=1.2，β0=0.020118.7103.914.8α0=1.2，β0=0.030115.5101.813.710.13245/j.hust.210921.F012图12不同阻尼系数下阻尼材料的优化分布4 结语研究了平稳激励下附加阻尼薄板结构的声辐射优化问题．建立有限元/边界元声振强耦合方程，推导并建立以阻尼层相对厚度为设计变量的人工阻尼材料模型，提出声功率和声压灵敏度的一种直接求解方法并研究了几种因素对声辐射优化的影响，得出了相应的优化结果．