声学覆盖层的有限元计算基本理论

1、第8章带圆柱空腔吸声覆盖层的有限元计算基本理论及有限元数值计算模型数值法与解析法的比较在吸声覆盖层声学机理的理论研究中，解析法较好地计算出覆盖层的吸声系数或隔声量（需要基体材料的实测参数），分析覆盖层的吸声或隔声机理，但是要完全用解析法去解决复杂结构问题，则是一件非常困难的事情。随着数值计算技术的发展功能强大的商用软件的出现，20世纪90年代开始，国外开始运用有限元法（TheFiniteElementMethod,FEM）分析具有双周期结构的Alberich声学覆盖层的平面波散射、反射特性，开发出声学覆盖层设计的有限元软件。在21世纪初，国内谭红波用有限元法对声学覆盖层的声特性进行了分析。王曼

2、对声学覆盖层用有限元方法进行了更深入的研究。有限元法可以比较好地解决非规则结构的声学性能计算问题，也可以研究声固耦合振动模态和振型，对吸声覆盖层的主要吸声机理有直观的认识。要设计出具有良好综合性能的吸声覆盖层，必须两者结合，以理论研究成果指导方向，用有限元法来具体设计。吸声覆盖层的特殊结构和特殊的应用环境，商用软件如ANSYS是不能直接用来计算的。基于Mead关于简谐波在周期性多耦合线性系统中的传播理论与Bloch周期性结构理论，是水中周期分布空腔结构体的有限元计算基本理论。1 Alberich吸声覆盖层的有限元数值计算模型20世纪90年代初，法国人Hennion和Decarpig印度人Eas

3、waran和Munjal发表了两篇关于Alberich吸声覆盖层的有限元数值计算模型的论文，提出了在用有限元方法计算共振腔吸声覆盖层时需要解决的几个共性问题，适用于任意形状空腔的谐振腔型吸声覆盖层。1 Alberich吸声覆盖层的有限元数值计算模型图1吸声覆盖层的单元结构1 Alberich吸声覆盖层的有限元数值计算模型 1)周期性结构处理 谐振腔型吸声覆盖层在x和y方向都无限延伸， z方向有空腔的存在而不均匀， 问题是三维的。 空腔沿x和y方向周期性或双周期性分布， 在作有限元分析时只取出其中的一个单元。1 Alberich吸声覆盖层的有限元数值计算模型 1)周期性结构处理 周期性用经典的B

4、loch关系式描述。 设x,y方向的周期分别为2d1和2d2， 设入射声压p， 与z轴的夹角为，1 Alberich吸声覆盖层的有限元数值计算模型 在xoy平面上的投影与x轴的夹角 。 Bloch关系式指出，简谐平面波( sincossinsincos)()-jk xyz iip x , y , z p e1 Alberich吸声覆盖层的有限元数值计算模型 简谐平面波作用下， 任何一个空间函数F（压力、位移等）都满足周期性关系式 此关系式在斜入射情况下特别有用。122sincos2sinsin12(2,2) = (, ) j d k-j d kF xd y d ,zF x,y z ee1 Al

5、berich吸声覆盖层的有限元数值计算模型 周期性结构的吸声覆盖层的反射系数和透射系数表示为 式中 和 分别是待确定的反射和透射系数；nmyjxjzjkPnmi ,eeeRx,y, zpx, y, zpmnnm,)()(nmyjxjzjkPnm,eeeTx, y, zpmnnm,)(PnmRPnmT1 Alberich吸声覆盖层的有限元数值计算模型 入射和透射波中包含有正常的传播波 又包含有z方向指数衰减的非传播波 。 高阶波是非传播波，可有限截断。122222/sincos ;/sinsin ;.nmnmnmndkmdkkk)0(2nmk)0(2nmk全反射全反射 当 透射波的方向偏向界面

6、临界入射角 当入射角 是虚数 声压反射系数的绝对值是1。全反射全反射 21cc112sincrcccr2sin1222cosnzz221sinsincc2121nnpnnzzRzz凋落波凋落波 全反射的折射波 expexpttyPAxi kyttxkiexpcosypAxkyt2expsinxyAvxkyt 2expcosyyyk Avxkyt 平均能量密度 平均能速度 能速度等于相速度，介于两种介质的声速之间 22222222222exp2exp242yyA kyAykc EykPvy非均匀波非均匀波 模式的波矢必须是实的 凋落波的波矢是复的 复波矢的平面波 expPAikr220 ckkR

7、Iikkk 两种情况 普通的平面波 非均匀波 和 垂直 ，凋落波 2202 RIRIRRIIRIiiickkkkkkkkkk220 RRIIckkkk0RIkk0Ik220 RRckk RkIk1 Alberich吸声覆盖层的有限元数值计算模型 2）包括流体耦合的吸声覆盖层单元的有限元处理 单元的有限元建模通常利用现有的各种有限元商用软件。 吸声覆盖层的一边是水， 另一边是钢板（单层或双层）。 吸声覆盖层的本体材料密度与水接近，水的耦合作用是重要的。1 Alberich吸声覆盖层的有限元数值计算模型 2）包括流体耦合的吸声覆盖层单元的有限元处理 对吸声覆盖层作有限元分析时必须包括一部分水单元。

8、设节点的压力和位移的向量是p和U， 单元的控制方程形式上可以写成 FpU222222cMLRcRMKsss1 Alberich吸声覆盖层的有限元数值计算模型 式中Ks和Ms是结构的刚度和质量矩阵； L和MP是水单元的压缩性和质量矩阵; R是反映结构和水耦合的耦合矩阵; F是作用力的节点值, 是流体域边界上的压力法向梯度的节点值. 水单元也可以用等效弹性单元来描述,取Lame常数。0,2c FpU222222cMLRcRMKsss1 Alberich吸声覆盖层的有限元数值计算模型 计算中包含的水单元究竟要多大值得研究. Easwaran和Munjal的工作说明， 水介质负荷主要在谐振频率附近起作

9、用, 偏离谐振频率后影响很小, 甚至可以不考虑水的作用。 这一点符合流固耦合的一般原理。 问题是在计算之初并不知道谐振频率, 合理大小的水单元必须包括在内.脉动球源 球状声源，作对称的周期的膨胀和收缩运动 体积变化速度，强度00exp()rrrVVit2004Qr V1011exprPikrcrikrV辐射的声场00000exp11VcrPikrrrikr000exp41iQikrrikrr000011exp11rV rikrVikrrrikr0011exp41Q ikikrikrrikrr 辐射的声强 总功率 222242100002222200022181PVrQfIcrcrkrcrkr

10、22242000022002121VrQfckrckr W 低频近似 ，波长比声源的尺度大 高频近似 ，波长比声源的尺度小 222200222832QfQckIcrr2202QfcW20222032QcIr r20208QcrW声场对声源的作用 声源表面受到的反作用00rrrFSP 00000exp()1ikrVcSitikr0rrFZ V 辐射阻抗 辐射阻抗 辐射阻 辐射抗0rrrrrrrFZRiXV 202202001rkrkcSRr2020001rkkrcSXr辐射对声源振动的影响 平均发射功率2012rRV 声源振动系统 考虑辐射对声源的作用 同振质量 0mFZ VmmKZRim0mr

11、FZZVrrXm 低频小球 ，辐射功率与频率平方成正比，辐射很弱 同振质量等于小球体积的介质的质量的3倍 高频大球 辐射抗和同振质量近似为零 2000rRcSkr000rXS r 00rRcS1 Alberich吸声覆盖层的有限元数值计算模型 计算中包含的水单元究竟要多大值得研究. Easwaran和Munjal的工作说明， 水介质负荷主要在谐振频率附近起作用, 偏离谐振频率后影响很小, 甚至可以不考虑水的作用。 这一点符合流固耦合的一般原理。 问题是在计算之初并不知道谐振频率, 合理大小的水单元必须包括在内.1 Alberich吸声覆盖层的有限元数值计算模型 3）边界条件的处理 单元四周和上

12、下两面的边界条件的处理是一个关键问题。 四周的边界条件取周期性条件， 在垂直入射情况下简化为平衡条件 与水接触的面的边界条件是压力连续 和振速（或压力的法向导数）连续。1 Alberich吸声覆盖层的有限元数值计算模型 3）边界条件的处理 在有限元计算中其它技巧。 利用1/4对称性大大缩短计算时间、 减少存储容量； 网格的划分可以按1/4横波波长作为准则等。1 Alberich吸声覆盖层的有限元数值计算模型 Alberich吸声覆盖层有限元数值计算结果 （1）计算结果与实验结果符合良好 用硅橡胶制作的Alberich吸声覆盖层的压力传输系数。 覆盖层厚度4cm，空腔高2cm， 直径1.5cm，

13、间距5cm。 材料的杨氏模量的实部E0=1.8106N/m2损耗因子E=0.15 泊松比v=0.49976，密度=1000kg/m3。1 Alberich吸声覆盖层的有限元数值计算模型 实线是实验结果， 虚线是Bennion的计算结果， 点线是Easwaran的计算结果。 两种有限元方法的计算结果稍有差别， 因边界条件的不同取法造成的。1 Alberich吸声覆盖层的有限元数值计算模型 谐振腔型吸声结构两种振动模式。 通过有限元计算描绘出谐振腔壁的两种谐振模式： 空腔壁的径向振动Alberich覆盖层在声波作用下的振动响应1 Alberich吸声覆盖层的有限元数值计算模型 谐振腔型吸声结构两种

14、振动模式。 表面覆盖层的似鼓形弯曲振动 有助于吸声机理分析和吸声结构设计。Alberich覆盖层在声波作用下的振动响应2水中周期分布空腔结构体的有限元计算基本理论吸声覆盖层敷设于水下目标壳体示意图2水中周期分布空腔结构体的有限元计算基本理论吸声覆盖层敷设于水中目标外表面，与目标外壳组成多层复合结构。(a）为单层壳体，整个多层复合结构分别由无限半空间水、声学覆盖层、耐压钢板、空气层组成。2水中周期分布空腔结构体的有限元计算基本理论双层壳体，则增加一个中间水层和外壳轻钢板层。吸声覆盖层厚度一般为厘米级，可以忽略外壳的弧度，认为是无限大的平面，整个结构可以简化为无限大的平面多层结构问题。2水中周期分

15、布空腔结构体的有限元计算基本理论吸声覆盖层中含有周期分布的空腔，有了“无限大”的假设后，取一个周期单元来进行分析，在平面波入射时，所分析的一个周期单元与其相邻单元之间的边界就可以利用Bloch周期边界条件来处理。2水中周期分布空腔结构体的有限元计算基本理论商用软件元法处理这种周期边界条件，可以利用ANSYS有限元软件划分单元节点，所有的计算均在Matlab环境下编程进行。2水中周期分布空腔结构体的有限元计算基本理论当声波在带有空腔结构的黏弹性材料中传播时，入射纵波在黏弹性材料与空腔的界面上发生剪切形变，必须将黏弹性材料作为具有一定内损耗的弹性介质来考虑，将黏弹性材料的杨氏模量和剪切模量表示为复

16、数，复弹性模量的虚部与实部之比，即为该模量的损耗因子。2水中周期分布空腔结构体的有限元计算基本理论关于边界条件。在声学覆盖层的无限大平面多层结构模型中，有多处边界条件：固体固体边界，如黏弹性材料与钢板之间，其边界条件是法向、切向应力分别连续，各方向质点位移分别连续；2水中周期分布空腔结构体的有限元计算基本理论固体流体边界，如黏弹性材料与水或钢板与水之间的界面，边界条件是固体中的法向应力与水中声压大小相等、符号相反，法向质点位移连续，切向应力和位移均为零；2水中周期分布空腔结构体的有限元计算基本理论固体空气边界，如黏弹性材料内的空腔表面或钢板后的界面，近似为自由边界，即法向和切向的应力均为零。2

17、水中周期分布空腔结构体的有限元计算基本理论声学覆盖层的声学性能计算是一个流固耦合力学问题。固体在流体载荷作用下会产生变形和运动，而这一变形和运动反过来又影响流场内的变化，改变流体中载荷的分布和大小。一般流固耦合方程的特点是方程的定义域同时有流体域和固体域，2水中周期分布空腔结构体的有限元计算基本理论未知变量既有流体变量又有固体变量，流体域和固体域通常无法单独求解，需要利用一定的显式或隐式方法求解方程组获取。2水中周期分布空腔结构体的有限元计算基本理论用有限元方法计算声学覆盖层声学性能的基本过程可简述如下。取一个多层周期单元，将入射端半无限水空间用一个平面S-（流体外边界面）划分成两部分，对与结

18、构相邻的这部分有限水层也取相应的一个单元对结构和流体区域分别建立运动方程和声场波动方程。2水中周期分布空腔结构体的有限元计算基本理论结构与流体边界处应力连续、位移连续，结构运动方程和声场波动方程是耦合方程。通过将结构和水区域分别进行有限元网格离散，形成结构和声场耦合的有限元方程。2水中周期分布空腔结构体的有限元计算基本理论在流体外边界面上，将声场作平面波展开，该面上的声压梯度转化为一个等效的节点载荷，作为流体有限元方程的载荷项。求解耦合的有限元方程得到流体节点声压，进一步求得各阶平面波的反射或透射系数。2水中周期分布空腔结构体的有限元计算基本理论 2.1流体有限元方程 声学覆盖层和钢板结构外是

19、半无限空间的流体水，满足均匀、不可压缩、无黏滞和流速为零的理想、流体条件，且声波激励是流体域内稳定状态附近的小扰动，则流体介质中波动方程矩阵形式表示为 式中c为流体介质中声速（ ，K为流体介质体积模量， 0为流体介质密度）;p为声压； 是梯度算子的转置22210Tppct0/Kc zyxT,2水中周期分布空腔结构体的有限元计算基本理论 取流体中一个六面体小单元。流体域中有声压的微小增量，对流体区域内积分p，得到 式中Vf是流体域；Ssf是流体、固体之间的交界面； Sff代表流体与流体之间的交界面，即将流体分隔为有限区域和无限区域的假想界面S-和S+; n表示表面的单位外法线向量。 2221ff

20、sfffTTVVSSppdVpp dVnpp dSct2水中周期分布空腔结构体的有限元计算基本理论 小振幅线性声场中介质运动方程为 式中u为流体介质的质点振速。 根据流固界面边界条件，设固体介质的质点位移为，在界面附近u=/t，流固界面附近固体质点位移与流体介质中声压的关系 写成矩阵形式为202pt 220tnpnTT0upt 2水中周期分布空腔结构体的有限元计算基本理论 式中 22011eeffesfTTppeeVVTpIeSNNdVpBB dVpccNnNdS eSeTpTIpeffdSpNnN2水中周期分布空腔结构体的有限元计算基本理论 固体结构振动的有限元方程比较常用，形式如下： 式中

21、Ks,Ms分别是固体结构的整体刚度矩阵、质量矩阵，分别由单元刚度矩阵、单元质量矩阵集成产生： eVTsseVTsesesdVNNKdVBDBMP = + sssMK FF2水中周期分布空腔结构体的有限元计算基本理论 结构有限元方程中右端的Fs项 是结构所受机械激励的等效节点载荷。 Fp是流体对结构作用的等效节点载荷。 由于流固界面上法向应力连续， 在界面上任一点流体对结构单位面积的作用力大小等于该点声压，方向沿界面的法线方向。2水中周期分布空腔结构体的有限元计算基本理论 力矢量fp的方向与流体区域在流固界面的外法向矢量方向一致，利用节点声压和面有限元的插值形函数可表示成 再将分布载荷fcp离散

22、作用到节点上，得到 IeTpIpnpNnpf 1113331311331333eeeeeeeNeTNpITNpINeTNpTNppTNppcpNnNnpNNfNf2水中周期分布空腔结构体的有限元计算基本理论 向量fep的元素是单元各节点的等效节点载荷各个方向的分量， 向量长度是结构节点的自由度总数。 对单元面积分再对流固界面上各单元求和， 就得到整个流固界面上流体对结构作用的总等效节点载荷，即 eSeTpIeSpepesfesfdSpNnNdSfF2水中周期分布空腔结构体的有限元计算基本理论 利用全体流体节点的声压向量p， 即 pdSNnNFeSTpIpesf pRFTp2水中周期分布空腔结构

23、体的有限元计算基本理论 最后得 sTssFpRKM2水中周期分布空腔结构体的有限元计算基本理论2.3 流体域外边界的边界条件?流体有限元方程中右端项的物理意义，是在人为划分的一个流体外边界面上，半无限流体域内声压作用在界面上的声压梯度的等效节点载荷。2水中周期分布空腔结构体的有限元计算基本理论一般考虑声学覆盖层后面也是无限流体情况，就存在声波的入射、反射和透射现象。先将流体的入射端面和透射端面的流体外边界面分别记为，两个面的z坐标值分别记为z-，z+,入射端面和透射端面声压梯度的等效节点载荷分别记为?-?+。2水中周期分布空腔结构体的有限元计算基本理论 由单元矩阵(H-)e集成得到总体矩阵H-

24、，则入射波产生的总节点载荷向量为 面S-上总的声压梯度的等效节点载荷向量为 S+面上透射波声压梯度的等效节点载荷为Hejkjkzicoscos pCiri2 pCN12水中周期分布空腔结构体的有限元计算基本理论 将流体质量矩阵MP、刚度矩阵KP、流体固体耦合矩阵R的元素所对应的流体节点的排列顺序与?中元素相对应的流体节点的排列顺序保持一致，整理得到 整个系统流固耦合有限元方程可以写为 0202 - - C = RpMCKpp 02022CFpRMCKMKRsppssT2水中周期分布空腔结构体的有限元计算基本理论 2.4 周期性边界条件 建立有限元模型时只取了个周期单胞，在该周期单元的边界上也必

25、须假定适当的边界条件。由于结构的周期性，在单元的S1与S2（或S3与S4）面上， 如果节点位置分布相同， 则两个面上对应流体节点的声压相同，仅差一个相位，忽略由于材料损耗而引起的两个面之间的位移幅度差，则两个面上对应的结构节点位移之间也有这样的关系。2水中周期分布空腔结构体的有限元计算基本理论2.5平面波反射系数与透射系数求解流固耦合有限元方程，可得到周期单元内部及边界上流体节点S1,S2,C1的声压及结构节点的位移。入射侧的流体节点反射波声压可由总声压减去入射波声压得到。2水中周期分布空腔结构体的有限元计算基本理论取流体外边界S-上N个节点，建立节点坐标与反射波声压关系的N个方程，求解方程即

26、得到各阶平面波的反射系数。如果透射端为水，则在透射端流体外边界上取N个节点，建立节点坐标与透射波声压关系的N个方程，求解方程得到各阶平面波的透射系数。2水中周期分布空腔结构体的有限元计算基本理论 在各阶反射波和透射波中，z方向波数knm为虚数的各阶波是在z方向呈指数衰减的非均匀波，对远场反射或透射系数没有贡献。 远场反射或透射系数仅由满足k2k2nx+k2my条件的各阶传播波决定。从能量观点出发，可用各阶波的反射或透射系数均方根求得远场反射系数或透射系数，即NTtNRRnmnmpnmnm/222水中周期分布空腔结构体的有限元计算基本理论 一般情况下只有（0,0）阶波为传播波的条件， 尤其在声波

27、正入射时，反射系数和透射系数就分别等于（0,0）阶波的反射系数和透射系数。NTtNRRnmnmpnmnm/223有限元方法编程实现与验证对于实际的吸声覆盖层结构，细管的穿孔系数一般很小，覆盖层底部喇叭空腔的穿孔系数最大可达0.5以上；同时计算的频率范围从几十赫兹到几十千赫兹。3有限元方法编程实现与验证若按空腔尺寸或声波波长来划分单元，则在顾及细管或高频时存储所需的内存太大，也影响计算速度，在有限元程序编制上可利用ANSYS划分单元，自编Matlab软件解决稀疏矩阵存储和计算，采用两者之长，可以基本解决这个问题。3有限元方法编程实现与验证 3.1Matlab中稀疏矩阵存储和计算 有限元方法是利用

28、计算机求解数学物理或工程实际中的微分方程问题或偏微分方程问题的一种系统化的数值计算方法。 有限元方法的一个显著特点是在迭代求解过程中涉及大量的矩阵计算， 可利用Matlab强大的数值和矩阵计算方面的能力和可视化图形功能来实现。 针对有限元编程中的稀疏矩阵和矩阵求解，Matlab提供了许多有用的函数。3有限元方法编程实现与验证 A=sparse(A)， 把总体系数阵A转化为稀疏阵； B=sparse(B)， 把总体系数阵B转化为稀疏阵； Result=AB.Matlab根据A情况，采用最优化方法 （Matlab提供了4种高斯消去法或其变种， 还提供了6个迭代法求解函数） 求解系数阵为稀疏矩阵的线

29、性方程组。3有限元方法编程实现与验证图5 Q矩阵非零元素位置示意图3有限元方法编程实现与验证 图是前面有限元模型计算中的Q矩阵内非零元素的分布情况， 这个矩阵是一个稀疏矩阵， 正好可以利用Matlab的各种矩阵计算函数。 如模型划分为1000个节点， 采用全矩阵存储和计算时，占用的内存约为4GB， 超过了个人计算机的内存极限，基本没有使用价值。3有限元方法编程实现与验证 而利用稀疏矩阵可以将占用内存降低到300MB。 在节点数继续增加的情况下， 内存节约的效果更明显。 在利用个人计算机进行吸声覆盖层模型计算时， 可以使模型节点数高达数万， 对穿孔率变化范围大的实际吸声覆盖层计算有着非常实用的意义。3有限元方法编程实现与验证3.2 有限元模型的单元类型选择使用ANSYS进行划分单元时所建的模型均为非对称的三维模型，几何拓扑结构复杂，且吸声结构几何形状复杂难以直接用六面体单元来划分网格，一般都采用适应性更强的四面体单元。对流体介质和吸声固体结构（橡胶，钢板等）采用三维结构单元（SOLID45）。3有限元方法编程实现与验证 它们都是三维八节点单元，需要先将它们转换为四节点单元。通过一些算例证实，