功能梯度耦合梁的能量辐射传递模型

王 幸， 钟 强， 李 翱， 陈海波

(中国科学技术大学 近代力学系 中国科学院材料力学行为和设计重点实验室，合肥 230026)

复合材料梁在航空航天领域具有十分重要而广泛的应用，如在飞机机翼、直升机尾桨等零部件上已大量使用[1]。传统层压复合材料结构在极端工况下，经常会出现脱层和裂纹等损伤失效现象。与传统的层压复合材料不同，功能梯度材料(functionally graded material, FGM)是一种新型的非均质复合材料，它是由两种或者两种以上的材料在微观上融合制备而成[2-3]。FGM中不同材料成分之间的过渡是渐变的，这可以有效降低分层的风险，所以传统层压复合材料梁越来越多的被FGM梁代替。像机翼、尾桨等大型部件中耦合梁结构是不可避免的，且在实际工况中往往受到高频载荷[4-5]，引发结构强度破坏和疲劳失效等问题。所以研究FGM耦合梁的高频振动是很有必要的。

发生高频振动时，子频带内的模态密集，模态高度重叠，且特征波长远小于结构尺寸。根据高频振动的特征，将带宽内模态数大于5或模态重叠因子大于2的频段定义为高频段。传统有限元(finite element method，FEM)[6]和边界元法(boundary element method，BEM)[7]在求解结构高频振动时会遇到两个问题：①由于振动波长较小导致的计算成本高；②结构不确定性非常敏感导致的计算鲁棒性不强。为了避免传统数值方法在高频段遇到的这两个问题，学者们提出了多种的能量分析方法。其中最流行的是统计能量分析法(statistical energy analysis，SEA)[8-9]。但受限于扩散场假设，SEA仅能给出振动子系统的平均能量响应，而Bot等[10-12]通过类比热辐射提出了能量辐射传递法(radiative energy transfer method，RETM)，它可以准确地预测系统中任意位置的能量响应。RETM理论遵循惠更斯原理与线性叠加原理，系统的振动场由域内实源和边界虚源叠加而成，每个源产生的波场对于一维系统是平面波，对于二维系统是柱面波，对于三维系统是球面波。但是目前RETM的应用局限于各向同性系统，本文将此方法拓展到复合材料梁的高频响应分析中。

在FGM梁的动力学研究方面，Aydogdu等[13]研究了简支边界条件下FGM梁的自由振动问题，赵亮等[14]针对轴向运动的FGM悬臂梁进行动力学分析，邓昊等[15]研究了沿轴向指数分布的FGM梁的固有频率和模态振型，Liu等[16]通过建立能量流模型计算了受横向激励作用下的FGM梁的高频振动响应。但是目前国内针对FGM耦合梁高频振动响应的研究较少。本文考虑FGM耦合梁中弯曲振动和轴向振动的耦合效应，基于RETM理论建立FGM梁振动响应的求解模型，通过数值算例，与波传播分析法[17](wave propagation analysis，WPA)计算的精确解对比，验证所建立的求解模型的准确性。WPA法是根据波的传播特性，通过对控制方程进行两次傅里叶变换建立的。第一次变换得到位移在波数空间的幅值，第二次逆变换得到位移在空间内的直接场解。位移解为直接场解加上一般解，一般解由边界条件确定，因此WPA方法相当于理论解，可以用于验证RETM的准确性。此外，本文还分析了材料梯度指数n、频率f以及阻尼损耗因子η对FGM梁频散关系及能量响应的影响，讨论了梯度指数n对FGM耦合梁能量传递系数的影响，分析了FGM耦合梁弯曲和轴向振动响应。

1 FGM梁的振动模型

1.1 材料特性

图1为矩形截面FGM梁，梁的长宽高为l×b×h，建立如图1所示的笛卡尔坐标系，材料在z方向按照幂函数形式分布

(1)

式中：P(z)为材料特性，包括弹性模量、质量密度和阻尼损耗因子等；Rp=Ptm/Pbm，Ptm和Pbm分别为FGM梁的顶层和底层的材料特性；n(n≥0)为材料梯度指数，表示材料组成成分的体积分布；当n<1时，顶层材料占主导；当n=1时，材料在厚度方向上的变化是线性的；当n>1时，底层材料占主导。

图1 矩形截面FGM梁示意图Fig.1 Schematic diagram of rectangular cross-section FGM beam

1.2 梁中的波传播

图2(a)为受横向点激励F(x,t)=δ(x-x0)F0ejωt作用的FGM梁，将产生横向振动。其中F0为激励的振幅，ω为圆频率，t为时间变量，j2=-1。参照Liu等的研究，基于Euler-Bernoulli梁理论，建立拉格朗日方程，再根据哈密顿原理，推导得到FGM梁的弯曲振动控制方程为

(2)

(3)

式中：w为梁的横向位移；A为梁的横截面面积；D1为拉伸刚度；D2为耦合刚度；D3为弯曲刚度；I1,I2和I3为惯性矩；E(z)为弹性模量；ρ(z)为质量密度。

图2 单个FGM梁的振动Fig.2 Vibration of single FGM beam

图2(b)为受轴向点激励F(x,t)作用的FGM梁，将发生轴向振动。不考虑弯曲振动对轴向振动的影响，可得轴向振动控制方程为

(4)

式中，u为梁的轴向位移。

假设行波解，w(x,t)=Anejkbxejωt，u(x,t)=Bnejklx·ejωt，kb和kl分别为弯曲波波数和纵波波数。分别代入式(2)和式(4)，解齐次方程得

(5)

横向位移和轴向位移的通解为

(6)

图3为无限FGM梁中的波传播示意图。图3(a)中弯曲波在FGM梁中有两个传播波和两个倏逝波，图3(b)中纵波包含两个传播波。倏逝波相对于空间呈指数形式衰减，传播了几个波长的距离后就可忽略不计。而在高频情况下，波长远小于结构尺寸，这意味着倏逝波仅在激励点附近存在。

图3 无限FGM梁中的波传播Fig.3 Wave propagation in an infinite FGM beam

考虑结构阻尼时，E(z)变为E(z)[1+jη(z)]，η(z)为迟滞阻尼损耗因子。有效阻尼ηeff用以描述在厚度方向上变化的阻尼对功率耗散的影响

(7)

1.3 能量密度与功率流强度

能量密度W=Wp+Wk，Wp为势能密度，Wk为动能密度。弯曲振动和轴向振动的时间平均能量密度〈W〉b和〈W〉l分别为

(8)

式中，〈·〉为时间平均；(·)*为取共轭；Re(·)为取实部。功率流强度为通过单位面积的功率。

弯曲振动和轴向振动的时间平均功率流强度〈I〉b和〈I〉l分别为

(9)

梁的时间平均输入功率〈Pin〉为

(10)

式中，Y为激励点的输入导纳

(11)

2 FGM梁的RETM模型

2.1 点源在自由场中的能量辐射

在使用RETM计算FGM梁振动响应之前，首先要明确RETM的有效域。开发RETM的初衷是用以解决阻尼较大时不满足扩散场假设导致无法用SEA计算的问题，所以理论上RETM的有效域是宽于SEA的。根据文献[18]，RETM相比SEA不需要对衰减系数和耦合强度进行限制，所以RETM的有效性条件为：子频带内的模态数N≫1，模态重叠因子M≫1。参照Bot等的研究，将N≥1和M≥1的共同部分视为RETM有效域。模态数用模态密度γ表示为N=γΔω，Δω为频段带宽；模态重叠因子定义为M=γηω，并且模态密度γ=L/(πcg)[19],cg=dω/dk为群速度。此外，根据式(5)不难推导得到：轴向振动的群速度等于相速度，弯曲振动的群速度与相速度正相关，相速度cp=ω/k。综上，相同频率下，波数越大，群速度越小，导致模态密度越大，进而模态数和模态重叠因子越大，即波数越大结构越容易满足RETM的使用条件。对于耦合梁，需要满足Ni≫1，Mi≫1，i=1,2。

RETM的基本假设：① 均匀线性振动系统，且系统处于稳态振动阶段；②迟滞小阻尼模型；③不考虑倏逝波引起的近场效应；④不考虑波传播过程中的相互干涉。

图4为FGM梁中无穷小单元的能量守恒。根据假设①，在稳态情况下，FGM梁中无穷小单元的输入功率等于输出功率与耗散功率之和

div·I+Pdiss=Pinj

(12)

式中：div·为散度算子；Pdiss为耗散功率；Pinj为输入功率。根据假设②，迟滞阻尼模型中的耗散功率和能量密度成正比

Pdiss=ηωW

(13)

图4 FGM梁中无穷小单元的能量守恒Fig.4 The energy conservation of infinitesimal element in FGM beam

将I=cgW和式(13)代入式(12)得

(14)

式中，m=ηω/cg为能量衰减系数，对于FGM梁，阻尼用式(7)中的有效阻尼。

令r=x-x0，式(14)的齐次解为

I=Ce-m|r|

(15)

所以，功率流强度为

(16)

将式(14)左右两边从负无穷到正无穷积分，图5描述了积分域，分为3个部分：(-∞,+ε)；(ε,+∞)；(-ε,ε)，ε→0(ε>0)其中前两部分积分为零，因为在这两个域内没有能量输入，即

(17)

求得C=Pin/2，进而功率流强度为

(18)

相应的能量密度为

(19)

2.2 单一FGM梁的RETM

根据假设④，忽略波传播过程中的相互干涉，域内点的总能量密度遵循线性叠加原理。总波场由内部实源产生的直接场和边界虚源产生的反射场叠加而成。图6为单个FGM梁RETM示意图，梁上任意一点M的能量由激励点处的实源ρS引发的直接场和两端点处的虚源σA和σB引发的反射场共同产生。结合式(18)和式(19)，M点的功率流强度和能量密度为

(20)

图6 单个梁的RETM示意图Fig.6 RETM schematic diagram of a single beam

式中，ρS=Pin/2。RETM只需求解出边界虚源强度就可以描述梁中任意点的功率流强度和能量密度，而边界虚源强度可根据边界处的功率流平衡条件求解。边界处的功率流平衡为

(21)

解得

(22)

(23)

将式(22)代入式(20)，梁上任意一点的功率流强度和能量密度即可求得。

2.3 FGM耦合梁的RETM

计算耦合梁的振动响应时，需要对耦合处的能量传递系数进行求解。图7为θ角度耦合的两个半无限梁，其振动包括弯曲振动和轴向振动的耦合，耦合处的能量传递系数为θ的函数。In为入射波，Re为反射波，Tr为透射波，下标b为弯曲波，l为纵波，e为倏逝波。梁1和梁2的位移为

(24)

根据耦合处的位移连续及力的平衡条件可解出式(24)中各个波的幅值[RebReb,eRelTrbTrb,eTrl]。规定当输入为弯曲波时，Inl=0；当输入为纵波时，Inb=0。各波的功率为

(25)

则耦合处的能量传递系数为

(26)

式中：Pc,i为梁i中c型入射波的功率；Pd,ij为梁i中入射波传递到梁j中产生的d型波的功率；τcd,ij为梁i中的c型入射波在梁j中产生的d型波的能量传递系数；当i=j时为能量反射系数，当i≠j时为能量透射系数，i,j=1, 2，c, d=b, l。

图7 θ角度耦合的两个半无限梁Fig.7 Two semi-infinite beams coupled with angle θ

图8为耦合梁的RETM示意图，在耦合点B处，分别包含梁1和梁2的边界虚源σB1和σB2。RETM理论中每个子结构内部的能量密度和能量流分布只由该子结构内部的实源和虚源决定，与其他子结构内部的虚源无关，且弯曲波和纵波相互独立。梁1和梁2中的任意点M的功率流强度和能量密度为

(27)

边界处的功率流平衡

(28)

图8 耦合角度为θ的耦合梁RETM示意图Fig.8 RETM schematic diagram of a coupled beams with a coupling angle of θ

3 数值算例分析

将WPA求解的能量响应作为解析解，用数值模拟验证RETM求解FGM梁的准确性。数值模拟中用到的材料组分如表1所示，顶层为钢，底层为氧化铝。FGM梁两端是简支的，在梁的中部施加简谐的点激励。首先讨论材料梯度指数n、频率f以及阻尼损耗因子η等变量对单个FGM梁振动响应的影响，然后考察耦合角度θ和n对FGM耦合梁能量传递系数的影响，最后分别计算FGM耦合梁弯曲振动和轴向振动的能量响应。在描述能量响应时，取能量密度的参考值为1×10-12J/m。

表1 FGM梁底层和顶层的材料特性

3.1 单一FGM梁

图2(a)所示的单个FGM简支梁，梁的长宽高：L×b×h=1.000 m×0.001 m×0.001 m，在x=x0处受到横向简谐点激励。

图9给出了不同梯度指数n下FGM梁弯曲波波数随频率f的变化情况。当n一定时，f越大，波数越多。因为波数定义为单位波长内完整波的数量，频率越大波长越短，因此波数越多。当f一定时，n越大，波数越少。因为n越大，底层材料氧化铝占比越大，弹性模量越大，则梁的刚度越大，从而波数越少。此外，由图9可以看出，高频情况下单一FGM梁弯曲波数的量级为1×102。传统FEA需要对结构进行单元离散，且要求每个波长至少用6个单元描述，这就需要大量的单元离散。因此，只需要两个自由度的RETM的小计算量优势就体现出来了。

图9 不同梯度梯度指数n的FGM梁的波数随频率的变化Fig.9 Variation of wave number with frequency of FGM beams with different gradient index n

图10给出了梯度指数n和η对输入功率Pm的影响，频率f=10 000 Hz。当f和n一定时，η越大，Pm越小，但是η的影响较小；当f和η一定时，n越大，Pin越大，且当n增大到一定值时，Pin几乎不随n的改变而改变。

图10 梯度指数n和η对输入功率Pin的影响，f=10 000 HzFig.10 The influence of gradient index n and η on the input power Pin, f=10 000 Hz

图11给出了梯度指数n和f对Pin的影响，阻尼ηim=ηbm=0.05。当η和n一定时，f越大，Pin越小。

图11 梯度指数n和f对输入功率Pin的影响,ηtm=ηbm=0.05Fig.11 The influence of gradient index n and f on the input power Pin, ηtm=ηbm=0.05

图12给出了梯度指数n对能量衰减系数m的影响，阻尼ηtm=ηbm=0.05，频率f=10 000 Hz。当f和η一定时，n越大，m越小，且当n增大到一定值时，Pin几乎不随n的改变而改变。

图12 梯度指数n对能量衰减系数m的影响，ηtm=ηbm=0.05Fig.12 The influence of gradient index n on the energy dissipation coefficient m, ηtm=ηbm=0.05

图13给出了梯度指数n为1的单一FGM梁弯曲振动情况下在频率和阻尼平面的RETM有效域，是由临界模态数线和临界模态重叠因子线围成的半无限有效域。而由图9，n越大波数越小，进而越难满足RETM的使用条件，即当n增大，有效域的临界模态数线将右移，临界模态重叠因子线将上移，反之亦然。

图13 单一FGM梁在频率和阻尼平面的有效域，n=1Fig.13 Validity domain of RETM for single FGM beam in frequency and damping plane，n=1

图14 单个FGM梁的能量响应水平Fig.14 Energy response level of a single FGM beam

图15给出了不同梯度指数n的FGM梁能量密度分布，阻尼和频率同图14的设置。由图15(a)可以观察到，所有情况下，RETM结果和WPA结果吻合都较好，从而验证了RETM求解FGM梁高频振动能量响应的准确性。图15(b)给出不同n下的RETM解，可以看出n越大的 FGM梁对应的能量衰减幅度越小，这符合图12的规律：当η和f一定时，n越大，m越小。此外，除了激励点附近，n越大能量密度越大，结合图10和11，n越大Pm越大，激励点附近由于衰减较少能量接近输入能量，而随着传播距离的增加阻尼的影响渐渐体现出来。

图15 不同n的FGM梁能量响应水平Fig.15 Energy response level of FGM beams with different n

图16给出了3组不同和下的FGM梁的能量响应，梯度指数n=1。由图16(a)和图16(b)知，当f一定，η越大，能量损耗越大，以致边界处的入射波和反射波幅值减小，波的传播路程减短，干涉波幅值减小，所以图16(b)中WPA解在激励点附近的振荡不如图16(a)明显。对比图16(b)和图16(c)，当η一定，f越大，能量损耗越大，图16(c)同样也有WPA解在激励点附近振荡幅度较小的规律。这是因为m=ηω/cg，η和f越大，m越大，进而边界处能量越小，WPA解的振荡幅度越小，RETM解和WPA解吻合得越好。

图16 不同η,f下FGM梁能量响应水平，n=1Fig.16 Energy response level of the FGM beam for different η and f, n=1

3.2 FGM耦合梁

考察图7所示的半无限FGM耦合梁，截面尺寸均为b×h=0.001 m×0.001 m，梁1、梁2的材料梯度指数相同且均为1，结构阻尼为0，频率为10 000 Hz。

图17为分别入射弯曲波和纵波时，能量传递系数随θ的变化。可以看出，总的能量传递系数恒为1，这符合能量守恒。对比图17(a)和图17(b)，两幅图分别有两条相对应的曲线是完全一样的，即τbl，11=τlb，11,τbl，12=τlb，12，这符合互易性。在图17(b)，当耦合角度接近90°时，τll，11达到最大，而τll，12接近于0。

图17 不同耦合角度的能量传递系数，n1=n2=1Fig.17 Energy transfer coefficients for various coupling angles, n1=n2=1

图18给出了弯曲波入射时的能量传递系数随和的变化：①τbb，11随着θ的增大而增大；②τbb，12随着θ的增大呈先减小后增大的趋势；③τbl，11随着θ的增大呈先增后减的趋势；④τbl,12随着θ的增大呈先增后减的趋势。n的影响主要集中在0～1，当n增大到一定程度时，能量传递系数几乎不随n的改变而改变。

图19为加载横向点激励的FGM耦合梁，梁1和梁2 的尺寸：Li×bi×hi=1.000 m×0.001 m×0.001 m，梯度指数：ni=5。在梁1的x0=0.5 m处受到横向简谐力，耦合角度θ=π/3。

图20为该耦合梁在频率和阻尼平面的RETM有效域，图20中临界模态数线和临界模态重叠因子线均只有一条线，这是因为梁1和梁2的材料、尺寸均相同。相比于弯曲振动的RETM有效域，轴向振动要求频率更高、阻尼更大，因为相同频率下，纵向波波数远小于弯曲波波数。耦合FGM梁的RETM有效域需要以轴向振动的有效域为标准。

图18 能量传递系数随n和θ的变化，f=10 000 HzFig.18 The effect of n and θ on energy transfer coefficients, f=10 000 Hz

图19 加载横向点激励的FGM耦合梁Fig.19 FGM coupled beams loaded by a lateral point force

图20 FGM耦合梁在频率和阻尼平面的RETM有效域Fig.20 Validity domain of RETM for FGM coupled beams in frequency and damping plane

图21分别给出了FGM耦合梁的弯曲振动和轴向振动能量响应，结构阻尼损耗因子：ηtmi=ηbmi=0.1，i=1,2，激励频率f=20 000 Hz，可以看出RETM解和WPA解析解吻合较好。由于梁1与梁2耦合角度不为零，在耦合处存在能量的跃变现象。图21(a)中耦合点附近弯曲振动能量的振荡只存在于梁1，而图21(b)，梁1、梁2的耦合处附近均存在明显振荡。这是由于弯曲波的波数较大，沿传播路径的能量衰减较大，导致梁2右端的反射波无法回传到耦合处，因此梁2的耦合处附近没有这种能量振荡。由于相同频率下，纵波的波数比弯曲波的波数少，沿传播路径的能量衰减较小，边界处的反射波能够传到耦合处，在梁1梁2的耦合处附近均存在明显振荡。

图21 FGM耦合梁的能量响应水平Fig.21 Energy response level of FGM coupled beams

4 结 论

引入RETM求解FGM梁的响应，并采用WPA计算的解作为精确解对比来验证RETM求解FGM梁高频响应的准确性。给出了FGM梁的RETM有效域，分析了梯度指数n、结构阻尼损耗因子η以及频率f对单个FGM梁能量密度的影响，讨论了梯度指数n和耦合角度θ对FGM耦合梁能量传递系数的影响，分别计算了FGM耦合梁的弯曲振动响应和轴向振动响应。计算结果表明：

(1) 梯度指数n越大，RETM有效性条件越难满足，表现为临界模态数线右移和临界模态重叠因子线上移，此外耦合FGM梁的RETM有效域要以轴向振动的有效域为标准。

(2)n越大，能量衰减幅度越小；η和f越大，能量衰减幅度越大。能量衰减幅度越大，WPA解的振荡幅度越小，RETM解与WPA解吻合越好。

(3) 总的能量传递系数恒为1，这符合能量守恒，且能量传递系数满足互易性。能量传递系数受θ和n共同影响，n的影响主要集中在0～1。

(4) 用RETM求解的FGM耦合梁的弯曲振动响应和轴向振动响应均可与WPA解吻合，说明了RETM可以用于求解FGM耦合梁模型，拓展了RETM的使用范围。并且发现在同一耦合系统中，相同频率下，纵波沿传播路径的能量衰减较弯曲波小。