杨明绥，王同庆，范真真

(1.北京航空航天大学 流体与声学工程实验室，北京 100083;2.中航工业沈阳发动机设计研究所，辽宁 沈阳 110015)

粘弹性材料具有良好的吸能、高阻尼特性，且与水能够产生很好的声阻抗匹配，被广泛应用于水下潜艇的吸声覆盖层和阻尼隔振层.因此，粘弹性结构的吸声、隔振特性研究，成为国内外学者的一个重要研究内容［1-7］.

在水下声散射数值计算与潜艇声隐身设计的研究中，大多采用声管测量得到的局域声阻抗模型进行描述，即认为吸声材料的表面声阻抗只与当地声质点法向速度有关，而与其他位置处声质点法向速度无关.然而，随着低频主动声呐的应用，声阻抗的非局域特性不容忽视，一些学者和研究机构已开展消声瓦大试样吸声特性的实验研究工作［3］.因此，如何建立非局域声阻抗模型，描述消声瓦声阻抗的非局域特性是十分必要的.在非局域声阻抗代数模型的研究方面，Faverjon等［8-9］进行了基础性的研究工作，其针对铝板－多孔介质－铝板多层系统的声透射问题，分别从实验和解析2个方面开展了等效声阻抗代数模型的研究.然而目前对于消声瓦等粘弹性结构，还没有关于非局域声阻抗特性及其模型的研究工作发表.

本文参考Faverjon关于透声的等效声阻抗模型，发展了能够适用于粘弹性材料表面声散射阻抗非局域特性描述的理论模型，采用粘弹性动力学有限元方法，开展吸声材料表面声散射阻抗的非局域特性研究.本文的研究将对对粘弹性吸声覆盖层优化敷设以及潜艇声隐身优化设计具有重要意义.

1 声阻抗矩阵模型

声散射阻抗是描述材料声吸收性质的一种有效方法，刚性背衬的粘弹性吸声材料声散射阻抗的物理模型如图1所示.设声散射表面的声阻抗为zn，表面法向振动速度为vn，入射声压为p.将粘弹性材料表面离散为i(i=1，2，…，N)个有限网格点.则在网格点i处的表面声压为pi，引起的表面法向振速及声阻抗值分别记为vni、zni.

图1 物理模型Fig.1 Physical model

1.1 局域声阻抗

局域声阻抗模型认为:当地的声阻抗值zlocn(ω)只与当地的声压p和表面法向质点速度vn有关，与其他位置处的声质点法向速度无关，表现了声阻抗的局域特性，可以用声管进行小试样测量.如下:

式中:ω为角频率.在每个网格点i处均满足方程(1)，即

则表面声压pi和表面法向振动速度vni满足:

1.2 非局域声阻抗

当声阻抗表现为非局域特性时，当地声压不仅与当地法向振动速度有关，还与其他点处的法向振动速度相关.此时，pi与vni满足:

1.3 声阻抗矩阵模型

注意到，方程(3)、(4)均可以写成如下通用形式:

式中:p为表面声压向量，vn为表面声质点速度向量，zn为声阻抗矩阵.由以上的分析可知，当声阻抗表现为局域特性时，声阻抗矩阵zn为对角阵;当声阻抗表现为非局域特性时，声阻抗矩阵zn为非对角矩阵，阻抗矩阵中阻抗值的大小、分布情况分别由材料的吸声性能、非局域特性决定，通常情况下zn为主对角占优矩阵.

同时，表面声阻抗矩阵zn与结构的传递函数矩阵H具有互逆性［8-9］，即满足

且

声载荷对粘弹性材料的作用，可以等效为表面离散点处的节点力.因此，当j点声压为pj，i点速度响应为vni时，代入方程(7)可求得传递函数Hij，进而可构建传递函数矩阵H，并最终得到声阻抗矩阵zn.

1.4 声阻抗代数模型

参照Faverjon关于声透射等效阻抗模型，本文针对刚性背衬下粘弹性板的声散射问题，建立了一个声散射阻抗模型:

式中:rij为i点与j点之间的表面距离;ζ为局域声阻抗值，且满足ζ=ζR+iζⅠ=znii=zlocn;ρR(rij，ω)、ρI(rij，ω)分别为声阻、声抗密度函数，满足:

由于刚性背衬和粘弹性材料声耗散的影响，声阻函数ρR(rij，ω)均为正值，本文采用余弦函数的绝对值进行描述.式中:LR、LI、λR、λI、φI是随ω变化的常数;LR、LI为控制幅值呈指数衰减的长度尺度参数; λR、λI为振动波长的控制参数;φI分别为声抗密度函数的相位控制参数.声阻抗值zij具有如下性质: 1)是频率ω的函数;2)与局域声阻抗zlocn(ω)密切相关;3)声阻抗值zij与表面距离rij有关，并且随着rij的增加，阻抗值zij呈指数衰减.当rij→0时，zij→zlocn;当rij→∞时，zij→0.并且有zij=zji.

可由下式求局域声阻抗值zlocn:

求相位控制参数φI:

然后由方程(8)、(9)曲线拟合得到参数LR、LI、λR、λI的具体数值和变化规律.

2 粘弹性有限元

本文采用粘弹性有限元方法，进行方程(7)中传递函数Hij的计算.粘弹性材料的复常数模量模型形式简单，在简谐激励下能较好地描述粘弹性材料的力学性能［4］.基于复常数模量模型，可以在频域内建立粘弹性(弹性)－声场的耦合振动方程［5］:

式中:M为粘弹性(弹性)材料的质量矩阵，MA为流体附加质量矩阵，CS、CA分别为结构阻尼矩阵、流体阻尼矩阵，δ、f分别为位移响应向量、外载荷向量，K为结构刚度矩阵.对于粘弹性材料，K为复矩阵;对于弹性材料，K为实矩阵.当流固耦合较弱时，忽略阻尼与流体附加质量的影响，方程(12)可转换为

式中，fm、fp分别代表结构外载荷向量、声载荷向量，质量矩阵M、刚度矩阵K满足:

式中:ρ为结构材料密度，Me、Ke分别为单元质量矩阵、单元刚度矩阵，Nδ为形函数矩阵，本文采用8节点等参单元.Bδ为应变矩阵，D为弹性矩阵，满足:

式中:v为泊松比.对于弹性结构，弹性模量E为实常数;对于粘弹性结构，E为复模量，满足

式中:ER、EI分别为复模量的实部、虚部;η为粘弹性材料的损耗因子，满足

复模量E与剪切模量G、体积模量K:

当入射声波已知，则散射表面各节点处的声压向量p均可知，此时可利用方程(13)求得表面各节点处的结构位移响应向量δ.由谐波假设条件可知，表面位移向量δ与表面速度向量v满足:

此时，表面法向振动速度向量vn与表面速度向量v满足:

式中:n为节点法向矢量矩阵.将vn代入方程(7)可求得传递函数矩阵H.

进行结构固有频率分析时，方程(13)可转化为

对于弹性结构，方程(21)转化为广义实特征值问题;对于粘弹性结构或弹性－粘弹性夹层结构，方程(21)为广义复特征值问题，此时φ是复特征向量，λ=ω2是复特征值.将复特征值λ表示如下:

则第n阶模态损耗因子为

式中:λR、λI是复特征值λ的实部和虚部.

本文采用Fortran自编有限元程序，计算刚性背衬上粘弹性材料的传递函数矩阵H，由方程(6)得到声阻抗矩阵zn.

3 算例

3.1 算例校核

为校核自编代码的数值计算精度.计算了文献［8］中弹性－粘弹性－弹性夹层板结构自由振动的各阶模态频率fn、模态损耗因子ηn，并与文献［8］中的结果进行了对比.矩形夹芯板的长宽尺寸为0.304 8 m×0.304 8 m，上、下弹性板的厚度为0.762×10－3m，粘弹性芯层厚度为0.254×10－3m.弹性板的密度为2.737×103kg/m3，弹性模量为6.89×1010N/m2，泊松比为0.3;粘弹性芯层的剪切模量为0.896×106N/m2，密度为0.999×103kg/m3，泊松比为0.5.

表1为本文结果与文献［8］结果的对比.较好的一致性证明:自编代码的数值精度基本能满足粘弹性材料结构传递函数矩阵H的计算要求.

表1 结果对比Table 1 Results comparison

3.2 数值模型

本文算例尺寸为500 mm×500 mm×10 mm，材料为3M ISD 112型粘弹性材料的粘弹性板，材料参数采用5参数ADF模型［9］进行曲线拟合得到.其中剪切模量G、损耗因子η的频域拟合值与实验测量值对比情况如图2所示.

图2 3M ISD 112材料属性Fig.2 3M ISD 112 material characteristics

刚性背衬位于z2表面，z1为声散射表面.计算中z2表面采用刚性边界，其余各面均为自由边界.模型简图如图3所示，粘弹性有限元网格如图4所示，表面离散点编号及分布情况如图5所示.

图3 模型简图Fig.3 Schematic illustration of the model

图4 粘弹性有限元网格Fig.4 Viscoelastic finite element mesh

图5 表面离散点的分布及编号情况Fig.5 Distribution of dispersing nodes and numbers on surface

3.3 声阻抗矩阵的计算

计算得到传递函数矩阵H后，可由方程(6)求得声阻抗矩阵zn.由于图5中1～40号节点为边界点，为忽略边界的影响，只计算了41～121号节点.图6给出了100、500、1 000 Hz时表面声阻抗矩阵zn的数值结果，其中行号、列号1～81，对应于图5中41～121号节点.

由图6可知，粘弹性板的表面声阻抗矩阵zn为对称、对角占优矩阵.由方程(4)可知，此时模型表面声阻抗具有较强的非局域特性，且声阻抗值zij与频率ω、表面距离rij均有关.

图6 100、500、1 000 Hz时，|zn|数值分布情况Fig.6 Distribution of numerical value|zn|at 100，500 and 1 000 Hz

3.4 非局域声阻抗的代数模型

由方程(10)可以求得表面局域声阻抗值zlocn.［100 Hz，1 000 Hz］时，局域声阻抗值zlocn的计算结果如图7所示.由图可知，局域声阻抗的幅值随频率升高而降低，并且局域声阻抗的声阻均为正值，声抗均为负值.

图7 局域声阻抗zlocn的计算结果Fig.7 Results of acoustic impedance

得到zn、z后，利用方程(8)可求得声阻密度函数ρR(rij，ω)、声抗密度函数ρI(rij，ω)随rij的变化曲线.然后采用非线性曲线拟合方法，得到方程(9)中各参数值.本文采用Origin软件完成非线性曲线拟合.表2为100、500、1 000 Hz的各拟合参数值.图8～10为各频率下声阻、声抗密度函数的计算值与拟合值的对比，由图可知曲线拟合具有较好的数值精度，说明了代数模型的合理性.为研究粘弹性材料表面声阻抗的非局域特性提供了有效的研究方法.同时由表2可知，随着频率的增大，波长控制参数λR、λI均减小，这与实际物理现象是一致的.且声阻、声抗密度函数随着距离的增加，均呈现出负幂指数的振荡衰减.

表2 拟合参数的数值Table 2 Value of fitted parameters

图8 100 Hz时，函数ρR，ρI的结果对比Fig.8 Results comparison of function ρR，ρIat 100 Hz

图9 500 Hz时，函数ρR，ρI的结果对比Fig.9 Results comparison of function ρR，ρIat 500 Hz

图10 1 000 Hz时，函数ρR、ρI的结果对比Fig.10 Results comparison of function ρR、ρIat 1 000 Hz

4 结束语

本文将Faverjon针对航空材料声透射问题所提出的阻抗模型，发展为适用于粘弹性材料声散射问题的非局域声阻抗模型，并将其应用于低频主动声呐探测时粘弹性材料表面声阻抗出现的非局域特性研究方面.

以粘弹性有限元为基础，针对粘弹性吸声材料的声散射问题，完成了声阻抗矩阵、非局域声阻抗代数模型的求解，分析了粘弹性板表面声阻抗的局域和非局域特性.结果表明:代数模型的计算值与拟合值具有很好的一致性.进一步说明了本文所提出的方法的正确性，且具有较好的计算精度.为研究声阻抗的非局域特性，提供了可行的数值计算方法.

同时，在下一步工作中将开展粘弹性材料声阻抗非局域特性的实验研究，进一步完成本文所提出理论模型和计算方法的验证.

［1］ROY P K，GANESAN N.A vibration and damping analysis of circular plates with constrained damping layer treatment［J］.Computers and Structures，1993，49(2):269-274.

［2］RAMESH T C，GANESAN N.Finite element analysis of cylindrical shells with a constrained viscoelsatic layer［J］.Journal of Sound and Vibration，1994，172(3):359-370.

［3］王曼.水声吸声覆盖层理论与实验研究［D］.哈尔滨:哈尔滨工程大学，2004:78-90.

WANG Man.Theoretical and experimental study on underwater anechoic coating［D］.Harbin:Harbin Engineering University，2004:78-90.

［4］李军强，刘宏昭，王忠民.线性粘弹性本构方程及其动力学应用研究综述［J］.振动与冲击，2005，24(2):116-121.

LI Junqiang，LIU Hongzhao，WANG Zhongmin.Review on the linear constitutive equation and its dynamics applications to viscoelastic materials［J］.Journal of Vibration and Shock，2005，24(2):116-121.

［5］邹元杰.水中阻尼复合壳体结构声振特性的数值分析［D］.大连:大连理工大学，2004:101-136.

ZOU Yuanjie.Numeric alanalysis of vibro-acoustic characteristics of underwater damped composite shells［D］.Dalian:Dalian University of Technology，2004:101-136.

［6］ZHANG Q J，SAINSBURY M G.The Galerkin element method applied to the vibration of rectangular damped sandwich plates［J］.Computers and Structures，2000，74(6):717-730.

［7］TRINDADE M A.Reduced-order finite element models of viscoelastically damped beams through internal variables projection［J］.Journal of Vibration and Acoustics，2006，128(4):501-508.

［8］FAVERJON B，SOIZE C.Equivalent acoustic impedance model.Part 1:experiments and semi-physical model［J］.Journal of Sound and Vibration，2004，276(3-5):571-592.

［9］FAVERJON B，SOIZE C.Equivalent acoustic impedance model Part 2:analytical approximation［J］.Journal of Sound and Vibration，2004，276(3-5):593-613.

［10］ZHANG Junjie，LI Tianyun，YE Wenbing，ZHU Xiang.Acoustic radiation of damped cylindrical shell with arbitrary thickness in the fluid field［J］.Journal of Marine Science and Application，2010，9(4):431-438.