粘彈性材料聲阻抗非局域特性的數值研究

楊明綏，王同慶，范真真

(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.