耦合板结构的非结构零阶能量有限元分析

周红卫， 陈海波， 王用岩

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



耦合板结构的非结构零阶能量有限元分析

周红卫， 陈海波， 王用岩

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

在模拟结构高频振动时，随着分析频率的增高，传统的数值分析方法，如有限元(FEA)和边界元(BEM)等，面临着计算成本急剧增大的瓶颈[1]。于是，以统计能量分析(Statistical Energy Analysis，SEA)和能量流分析法(Energy Flow Analysis，EFA)等为代表的高频分析方法已经越来越广泛地应用于结构的高频振动及噪声分析中。

目前广泛使用的SEA基于模态理论，主要研究变量是耦合结构子系统的总能量；而EFA基于波动理论，其研究的主变量为波长和周期上双重平均的能量密度，在中高频内能较为准确地反映能量响应水平在空间上的分布，相比于SEA能得到更为丰富的结果，具有一定的优越性。Nefske等[2]利用有限元技术求解了梁内弯曲波的EFA偏微分方程，这种求解格式被称为能量有限元(Energy Finite Element Method，EFEM)。后续不少学者将EFEM扩展到杆、梁[3-4]、板[5-6]等常见基本结构中。2000年Wang[7]提出了零阶能量有限元法(EFEM0)，这种方法结合了SEA和EFEM的优点[8]，它利用有限体积法离散EFA偏微分方程，最终可以得到与SEA类似的线性方程组。但是，迄今EFEM0均基于结构网格，需要在板结构上划分结构化的四边形网格，这在外形不规则时往往会出现畸形单元，影响计算精度。如以三角形划分几何外形，即可克服上述问题。因此，发展非结构网格的零阶能量有限元法(Unstructured Zero-Order EFEM，uEFEM0)显得非常有实际意义。

基于以上原因，本文推导了三角形网格的uEFEM0的离散格式。对于耦合结构，结合考虑多波场时计算耦合板结构连接处耦合矩阵的一般处理方法[9]，对L型耦合板及简化汽车外壳的弯曲、拉伸和剪切波场的能量响应进行了数值模拟，并利用EFEM及商业软件AutoSEA对相同结构分别进行了分析，与本文计算结果进行了对比验证。

1非结构零阶能量有限元法

1.1薄板中的能量流分析方程

在高频振动的薄板中，存在三种不同形式的行波：弯曲波、拉伸波和剪切波。不同频率和不同种类的行波都互不相关地传播，并且它们的能量密度分布满足以下方程[10]：

(1)

式中:e是空间和时间上双重平均的能量密度，πin是单位面积上的输入能量，η是阻尼损耗因子，cg是板中行波的群速度，ω为圆频率。下标B、L和S分别代表弯曲波、拉伸波和剪切波。

1.2非结构网格有限体积法离散

刘学哲等[11]推导了辐射扩散方程的非结构有限体积法格式。本文基于此，推导了式(1)的离散格式。式(1)中的三种行波的控制方程形式一致，下列推导过程中略去了下标B、L或S。

图1　相邻控制体Fig.1 Adjacent control volume

将控制方程(1)在控制体i上进行积分：

(2)

由Gauss散度定理，有：

(3)

(4)

(5)

(6)

图2　共同边界上定义的向量Fig.2 Vectors on the common boundary

另外，认为e在Vi内恒定，则有：

∫ViηωedV=ηωeVi

(7)

并在控制体i内积分输入能量为：

∫Viπin,idV=Πin,i

(8)

综合式(5)、(7)和(8)，可得uEFEM0在控制体i上的单元方程：

(9)

其中:

(10)

可以验证，如果划分规则的四边形网格，由uEFEM0得出的计算格式(9)和Wang[7]推导出的应用结构化网格的EFEM0的计算格式是完全一致的。

1.3耦合边界上的能量流分析

Bitsie[12]推导出在无能量损耗的耦合边界上，能量流与各板在耦合处的能量密度之间的关系如下

(11)

式中:qcp为耦合边界上的能量流向量，I为单位矩阵，τ为功率传输系数矩阵，cg为由群速度组成的对角阵，ecp为各板在耦合边界上的能量密度。

若考虑编号为i和j的两块板耦合，并考虑三种行波的情形，则有如下具体形式：

(12)

(13)

由于不同子系统耦合节点上的能量密度是不相等的，因而在有限元离散的EFEM模型中需要增加耦合边界上的额外节点，称为重节点。而在EFEM0中，能量密度被设置在单元中心处，耦合边界上并无节点，因此无须设置重节点，而只需将qj在耦合边界上的计算格式相应地从式(5)改为式(11)，而后组装到总体刚度方程即可。由于这一原因，只需找出FEM模型中的耦合边界，便可将此模型直接用于EFEM0的后续分析。

1.4总体刚度方程的具体形式

首先将式(9)表示的单元刚度方程分别代入三种行波的控制方程式(1)，可以得出不含耦合信息的系统总体方程：

(14)

式中:

eP={e1Pe2P…eNP}，P=B,L式S

(15)

N为系统的控制体个数，BB等为相应波形的总体刚度矩阵。

接着找出系统中所有的耦合单元对，并应用式(11)添加耦合边界处各波场能量相互转化对于总体刚度矩阵的影响，得到最终的总体刚度矩阵为：

(16)

2算例

2.1L型耦合板

L型耦合板(如图3所示)的材料均为铝：密度为2 700 kg/m3，弹性模量为71 GPa，泊松比为0.3。假设三种波场全分析频率范围内的内损耗因子均为0.05。两板尺寸均为1 m×1 m，板1厚度为5 mm，板2厚度为1 mm。在板1中心处加载1 W的弯曲波能量，加载频率范围为1 kHz～10 kHz的11个1/3倍频程，各1/3倍频程中心频率fc如表1所示。用商业软件ANSYS划分了如图4所示的三套网格，称为Mesh1～Mesh3，各模型的单元数分别为100、400和444。

表1　各1/3倍频程的中心频率

图3　弯曲波输入的L型耦合板Fig.3 L-shape plates with bending wave load

图4　L型板的三套网格Fig.4 Three mesh models of L-shape plates

图5显示了Mesh3在1/3倍频程No.9(fc=6.3 kHz)内各波场的能量密度分布。各波场能量密度在两板的耦合边界处均不连续。对于弯曲波能量密度，加载点处最高，距加载点越远则越小；而面内波能量密度，耦合边界处最高，板2中面内波能量比板1大。这是由于该结构的外界能量输入只有板1中心的弯曲波能量，而面内波是来自于弯曲波传播到耦合边界上时的反射和折射,且面内波对弯曲波的透射系数远大于反射系数。

图5　耦合板各波场的能量响应；参考值：1x10-12 J/m2Fig.5 Energy density of L-shape plates; ref=1x10-12 J/m2

图6中为三种网格下的相同位置(如图3中的虚线Line1,1/3倍频程No.11)各波场能量密度的对比情况，横坐标为Line1上距点A的距离，点A、点B和点C为各边的中点。可以看出，uEFEM0具有良好的网格收敛性。另外，在乱序网格Mesh3上良好的适应性也说明了uEFEM0能适用于不规则结构的三角形网格划分。Park等[13]通过算例指出，随着分析频率的增高和内损耗因子的增大，在板2末端的面内波能量密度可能超过弯曲波能量密度，图6取得了与之类似的结果。这就使得在阻尼偏大或频率较高的情况下，考察距离加载点较远的振动结构能量响应时，面内波往往不能忽略。

图6　三套网格的相同位置的能量密度结果对比Fig.6 Comparison of energy density result of three mesh models on the same position

另外，对该结构用常规EFEM(划分如Mesh2网格相同尺寸的四边形网格，仅计算了弯曲波)以及商业软件AutoSEA分别进行了分析。三种方法得到的各板总能量结果对比情况(Mesh2)如图7～9所示。从各图中可以看出，EFEM和uEFEM0的弯曲波计算结果基本完全吻合，这说明本文提出的uEFEM0计算精度与传统EFEM一致。另外，uEFEM0和AutoSEA得到的各波场结果吻合较好。板2的弯曲波能量有1.5 dB左右的差距，板1的高频拉伸波结果有1～1.5dB左右的差距。这是由于能量有限元和统计能量法两种方法的差异造成的。相比较游进等得到的结果[6]，本文考虑了全部三种波形，拉伸波的对比结果更好一些。

图7 L型耦合板弯曲波总能量uEFEM0、EFEM和SEA结果对比Fig.7ComparisonoftotalbendingenergyresultofL-shapeplates,byuEFEM0,EFEMandSEA图8 L型耦合板拉伸波总能量uEFEM0和SEA结果对比Fig.8ComparisonoftotallongitudinalenergyresultofL-shapeplates,byuEFEM0andSEA图9 L型耦合板剪切波总能量uEFEM0和SEA结果对比Fig.9ComparisonoftotalshearenergyresultofL-shapeplates,byuEFEM0andSEA

通过对比，能说明本文处理的正确性。另外，相比于SEA，能量有限元能得到如图5～6所示的空间分布的能量结果，有助于工程上找到结构的危险位置。

2.2简化的汽车外壳

图10　简化的汽车外壳的外形及其uEFEM0模型Fig.10 uEFEM0 model of simplified vehicle shell

如图10所示，简化的汽车外壳材料为铝(参数同L型耦合板算例)，各板厚度均为5 mm。uEFEM0模型的网格数为468，结构中共有耦合单元79对。分析的频率范围为1 kHz～10 kHz(如表1)，假设全频域内的各波场阻尼损耗因子均为0.06。车体尺寸为：车身长度3 m，宽度1.5 m，高度1 m；车头高度0.4 m，顶盖长度2 m。设4个车轮在底板坐标系XOY内的坐标为(0.6,0.3)、(0.6,1.2)、(2.5,0.3)和(2.5,1.2)。该汽车壳侧板形状不规则，用结构网格不易划分,但易于划分非结构化网格。该案例可以初步体现非结构网格在处理不规则形状时的优势。

在车身底板4个车轮处施加弯曲波能量，假设分析频域内的4个外部功率输入均为0.25 W。图11从两个角度显示了结构在中心频率为fc=6.3 kHz的1/3倍频程内的弯曲波能量密度分布。图中的单位为dB，参考值为1x10-12J/m2。可以看出，底板上加载激励的车轮处的能量最高，而据底板最远的顶板处能量最低。

图11　简易汽车外壳的弯曲波能量密度分布；fc =6.3 kHzFig.11 Bending energy density distribution of simplified vehicle shell. fc =6.3 kHz

用AutoSEA建立对应的SEA模型，并在底板子系统全分析频段内加载弯曲波能量1 W。图12为全分析频段内汽车底板的总能量结果对比。从对比结果中可以看出，两种方法的弯曲波能量结果很接近，面内波能量有大约1～2 dB的差距，这和L型耦合板算例内取得的结果类似。

图12　底板的各波场总能量uEFEM0及SEA结果对比Fig.12 Comparison of total energy result of bottom plate, by uEFEM0 and SEA

3结论

本文基于薄板结构中三种波场EFA控制方程，推导了应用三角形网格的非结构网格零阶能量有限元的离散格式，提出了用于求解薄板高频振动的非结构零阶能量有限元(uEFEM0)分析方法。并运用uEFEM0计算了两个算例，均取得了与EFEM和AutoSEA较为一致的结果，验证了本文提出的方法的正确性。零阶能量有限元在耦合边界上不需要设置重节点，简化了使用FEM模型进行高频能量流分析的步骤。相比于EFEM0，由于uEFEM0采用三角形网格，能适应更为复杂的结构外形，具有更为广阔的工程应用前景。

另外考虑到三维声空间易于用四面体网格划分，使得本文方法在考虑声振耦合方面具有较大的优势。

参 考 文 献

[1]伍先俊, 朱石坚, 曹建华. 结构声振研究的功率流方法[J]. 力学进展, 2006, 36(3): 363-372.

WU Xian-jun, ZHU Shi-jian, CAO Jian-hua. Review on structural vibration power flow prediction methods[J]. Advances in Mechanics, 2006, 36(3): 363-372.

[2]Nefske D J, Sung S H. Power flow finite element analysis of dynamic systems: basic theory and application to beams[J]. Journal of Vibration Acoustics Stress and Reliability in Design, 1989, 111(1): 94-100.

[3]Wohlever J C, Bernhard R J. Mechanical energy flow models of rods and beams[J]. Journal of Sound and Vibration, 1992, 153(1): 1-19.

[4]孙丽萍, 聂武. 杆, 梁结构振动能量密度的简化解法[J]. 哈尔滨工程大学学报, 2004, 25(4): 403-406.

SUN Li-ping, NIE Wu. Asimplified method for solving energy density of rods and beams[J]. Journal of Harbin Engineering University, 2004, 25(4):403-406.

[5]Bouthier O M, Bernhard R J. Simple models of the energetics of transversely vibrating plates[J]. Journal of Sound and Vibration, 1995, 182(1): 149-164.

[6]游进, 李鸿光, 孟光. 耦合板结构随机能量有限元分析[J]. 振动与冲击, 2009, 28(11): 43-46.

YOU Jin, LI Hong-guang, MENG Guang. Random energy finite element analysis of coupled plate structures[J]. Journal of Vibration and Shock, 2009, 28(11):43-46.

[7]Wang Shuo. High frequency energy flow analysis methods: Numerical implementation, applications, and verification[D]. Purdue University: United States, 2000.

[8]Moens I. On the use and the validity of the energy finite element method for high frequency vibrations[D]. Catholic University of Louvain,2001.

[9]Langley R S, Heron K H. Elastic wave transmission through plate/beam junctions[J]. Journal of Sound and Vibration, 1990, 143(2): 241-253.

[10] Dong J, Choi K K, Wang A, et al. Parametric design sensitivity analysis of high‐frequency structural-acoustic problems using energy finite element method[J]. International Journal for Numerical Methods in Engineering, 2005, 62(1): 83-121.

[11] 刘学哲, 余云龙, 王瑞利, 等. 非结构任意多边形网格辐射扩散方程有限体积格式[J]. 数值计算与计算机应用, 2010,31(4): 259-270.

LIU Xue-zhe, YU Yun-long, WANG Rui-li, et al. A cell-centered finite volume scheme for discretizing diffusion equation on unstructured arbitrary polygonal meshes[J], Journal on Numerical Methods and Computer Applications, 2010, 31(4): 259-270.

[12] Bitsie F. The structural-acoustic energy finite-element method and energy boundary-element method[D]. Purdue University: United States, 1996.

第一作者 周红卫 男，博士生，1987年6月生

摘要：基于薄板高频振动时三种波场的能量流控制方程，用三角形网格划分板结构，推导了非结构零阶能量有限元的计算格式。并运用该方法计算了L型耦合板和简化的汽车外壳，得到各波场在空间上的能量密度结果。将计算结果与统计能量法和能量有限元进行对比,验证了其正确性。计算结果表明，频率很高时面内波能量往往能达到弯曲波能量水平，高频振动仿真时有必要同时考虑三种波场。

关键词：非结构零阶能量有限元；耦合板结构；多波场高频振动

Unstructured zero-order energy finite element method for coupled plate structures

ZHOUHong-wei,CHENHai-bo,WANGYong-yan(CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department of Modern Mechanics, University of Science and Technology of China, Hefei 230026, China)

Abstract:Based on the governing equations of energy flow analysis (EFA), meshing plates with triangles, an unstructured zero-order energy finite element method (uEFEM0) was developed. A procedure used to calculate L-shape coupled plates and a simplified vehicle shell was presented considering both bending and in-plane wave fields. The proposed method was used to predict the distribution of energy response of plates. To confirm its validity, the energy finite element method (EFEM) and statistical energy analysis (SEA) were employed to simulate the same structures, and the results agreed well with those using the proposed method. It was shown that in order to simulate plates vibration within a higher frequency range, it is necessary to consider not only the bending wave field but also the in-plane wave field, since the in-plane wave energy level may be close to the bending wave energy level.

Key words:unstructured zero-order finite element; coupled plates structure; high frequency vibration of multi-wave fields

中图分类号:TU311.3

文献标志码：A

DOI:10.13465/j.cnki.jvs.2015.13.024

通信作者陈海波 男，教授，1968年2月生

收稿日期：2014-03-03修改稿收到日期：2014-07-23