基于有限元的功率流分析方法及实现_谢基榕

第13卷第1期船舶力学Vol.13No.1 2009年2月Journal of Ship Mechanics Feb.2009文章编号：1007-7294（2009）01-0144-06
基于有限元的功率流分析方法及实现
谢基榕，吴文伟，沈顺根
（中国船舶科学研究中心，江苏无锡214082）
摘要：功率流分析在结构减振降噪研究中具有重要意义，但目前用于分析复杂结构功率流的工具非常缺乏。

鉴于有限元方法目前已广泛应用于结构动力学分析并且有许多成熟的商用软件，本文在Msc/Nastran的基础上二次开发了结构功率流分析模块并将其整合到Nastran的系统求解器中，使之成为工程人员分析结构功率流的专业工具。

关键词：功率流；有限元；DMAP
中图分类号：TB53文献标识码：A
Power flow analysis based on FEM:
theory and implementation
XIE Ji-rong,WU Wen-wei,SHEN Shun-gen
(China Ship Scientific Research Center,Wuxi214082,China)
Abstract：Analyzing power flow transmission in structures is of significance to control vibration and ac-companing noise,but present tools to analyze power flow in complicated structures is very lacking.This pa-per is concerned in redeveloping module for structure power flow analysis based on Msc/Nastran and inte-grating it into Nastran’s solve sequence system.This module will be a powerful tool for engineer to analyze power flow.
Key words:power flow;FEM;DMAP
1前言
机械振动引起结构辐射噪声是工业生产的一大污染源，于是减少振动并降低由此引起的辐射噪声便成了一项很重要的研究内容。

有限元法是目前应用于结构分析方面很重要很成熟的方法之一，它适合于分析结构静力及低频段的结构动力问题，目前已被广泛应用于工程研究领域。

结构振动传递，本质上是结构中振动能量的传播，如果能准确描述振动能量在结构中的传播路径并采取合理措施，那么在治理振动噪声上便能收到事半功倍的效果。

目前有限元方法已被广泛应用于结构动力学分析领域，并形成了许多成熟的商业软件，如Nastran，Ansys，Abacus等，这些商用软件通用性和二次开发功能很强，可在其基础上自主开发具有功率流分析功能的程序模块。

伍先俊[1]在使用有限元方法进行功率流分析方面进行过研究，他应用Ansys进行结构振动分析并将Ansys形成的结构有限元单元矩阵数据及位移响应解导出，然后使用Matlab将结构的单元矩阵组装成结构的广义刚度矩阵并结合位移计算出结构的作用内力，进而可计算出结构任意界面上的功率流。

然而，依照该方法分析结构功率流需要耗费大量人力，因为从Ansys中导出有限元矩阵数据，到计
收稿日期：2008-08-18
作者简介：谢基榕（1978-），男，中国船舶科学研究中心博士研究生，高级工程师。

算结构功率流需要用到不同软件，尤其是在功率流分析界面划分时，需要人工根据界面进行单元分组处理，不便于工程应用。

本文应用有限元分析软件Msc/Nastran 提供的开发语言DMAP 将功率流分析功能整合到Nastran 的系统求解器中，便可应用常规的Nastran 使用方式进行结构功率流分析，并在Patran 中处理计算结果。

2功率流计算方法
结构有限元的动力方程在激励为单频稳态力时可以写成
-ω211M +j ω11B +11
K 11ωωu =ωωF (1)
式中，11M ，11B ，11K 分别为结构的质量矩阵、阻尼矩阵及刚度矩阵，ωωu ，ωωF 分别为结构频率响应的位移解向量和激励力向量，位移解及激励力均使用复幅值表示。

根据功率的基本定义，计算结构某一界面处的功率流动需要先确定结构在该界面处的速度响应及界面处的内力。

事实上，在完成对整体模型的频率响应分析后，就已经获得了结构所有有限元节点处的速度响应，但在有限元分析软件中一般不输出单元间的内力数据
，因此要求根据结构单元特性数据及响应解计算出内力，才能进一步计算结构界面间的功率流。

从方程（1）式中可以看出，由结构单元的质量矩阵、阻尼矩阵及刚度矩阵组成的广义刚度矩阵
（R =-ω2
11M +j ω11B +11
K 11）联系了模型节点的位移解与作用于节点的外激励力之间的关系。

如果将模型单元根据待分析功率流的界面将模型单元划分成两组,如图1所示—待求功率流的界面S 把
模型单元G 分为两组（G 1∪G 軍1），则此时G 1
组模型单元的广义刚度矩阵关联了位移响应解与作用于
G 1节点的“外力”，而此“外力”则由分析界面的内
力与实际作用于该组单元的外力组成。

因此，可以计算出通过界面S 各节点流入单元组G 1的功率流。

假设全模型在外力系作用下的稳态响应解为
ωωu ，G 1组单元组成的广义刚度矩阵为R 11
1（仍将全模型的节点都考虑在内，则该刚度矩阵的阶数与全模型的响应阶数一致），从而可以得到在S 界面上作用于G 1组单元的力向量。

F 1ωω=R 1ωωωωu
(2)
将此力向量与速度响应向量j ωωωu 11
的共轭向量逐个元素相乘ωωP i =F 1ωωi
-j ωωωu i 11，（i =1,2,3,…,6N ）
(3)
即可得到流经第i 个自由度进入G 1组单元的功率流。

可以肯定，在ωωP 向量中只有S 界面上节点及G 1组上实际外力作用节点的自由度对应的元素才取非零值。

根据上述分析，可以任意取单元分组界面进而计算得流经任意自由度进入任意单元的功率流。

如果根据单元连接关系将全模型单元划分成取若干个单元分组使得在这些分组内任意两个单元都没有公共节点，那么可以求得经过任意节点流入任何单元的功率流。

从理论上说，无关分组数不大于节点关联的最大单元数。

图1界面S 将模型单元分成G 1与G 軍1两组Fig.1Interface S divides model into two groups
G 1
and G
軍1
第1期谢基榕等：基于有限元的功率流分析方法 (145)
146船舶力学第13卷第1期
3Msc/Nastran求解序列
Msc/Nastran软件将其求解代码分为两个层次[2]，首先，由基本计算机编程语言开发成实现一定目的的功能模块，由这些功能模块组成Nastran底层调用库，在功能模块库的基础上，Nastran用DMAP 语言编写求解序列调用模块库提供的基本功能模块，从而完成复杂问题的结构化求解过程。

DMAP语言本身也是一种结构化语言，可以用多个DMAP指令序列构成具有一定功能的子程序（SubDMAP），DMAP指令提供了参数操作、运行控制及模块（或子程序）调用功能。

在Nastran系统中，调用功能模块库完成有限元分析过程的完整的DMAP指令序列被称为求解器。

Msc/Nastran软件在发行时，系统自带通用问题的求解器及其DMAP源码，并且提供了用户改写求解代码的机制，从而为用户提供了求解功能扩展的基础。

Nastran用于进行结构频率响应直接法分析的是108号求解器，它的DMAP主程序叫SEDFREQ，由其调用更多的子程序完成求解过程中不同阶段的功能。

调用过程为：
誗“SEDFREQ”，结构频率响应分析直接法求解序列的主程序，本块代码对不同SubCase等条件进行循环控制，循环内容包括有限单元的数据准备阶段、方程求解及数据恢复等。

—“SUPER1”，有限单元的数据准备，包括读入数据卡片、单元矩阵形成、矩阵组装消减等。

*“IFPL”，读入数据文件，并将数据卡片内容排序组成内部数据块。

*“PHASE0”，对输入数据块进行整理操作。

*“PHASE1DR”，完成有限单元的数据准备，形成‘G’集矩阵并进一步消减到‘J’集。

·“PHASE1A”，有限单元的单元数据准备，形成‘G’集总体矩阵。

·“……”，其他一些矩阵消减操作。

—“FREQRS”，完成频率响应分析，是直接法频率响应分析的核心子程序，该子程序根据结构特性是否频率非线性进行分别处理。

当结构特性具有频率非线性时，对待分析的频率逐个生成结构单元矩阵并组装消减，在完成各个频率的响应解后再将所有解合并成单个数据块。

—“VDR1”，处理XYPLOT数据请求。

—“DISPRS”，在不同自由度集中恢复位移解，含有‘P’集位移解，该集包含‘G’集自由度。

—“SUPER3”，恢复其他解数据。

从108号求解器的求解序列代码布局上看，为了获得各节点的位移及单元矩阵数据，与之相关的子程序是“PHASE1DR”、“DISPRS”，当结构特性具有频率非线性时，“FREQRS”子程序也生成辅助的单元矩阵数据，从而这三个子程序将是我们研究的重点内容。

4修改Msc/Nastran求解器
根据第2节的理论分析可以很容易地确定分析功率流的思路，即首先必须准备全模型各节点的位移解，然后分别取不同的部分单元组，组装相应单元的质量矩阵、阻尼矩阵及刚度矩阵，配合全模型的位移解计算得作用该组单元的“内力”，进而计算得流入该组单元的功率流。

从程序实现的角度看，完整的功率流分析过程可以分为全模型位移求解及功率流回代两个阶段。

首先，对全模型单元进行频率响应分析，但是在计算得各节点‘G’集位移解后将此数据保存到外部文件中，完成位移解数据准备阶段；然后，取部分单元组进行频率响应分析，当求解进行到完成‘G’集矩阵时，将‘G’集矩阵保存到临时的外部文件后返回正常的频率响应分析求解序列；稍后，在本应该计算位移响应解的时候将保存于外部文件中的‘G’集的位移数据及矩阵数据读入，计算作用于单元“内力”及流经各节点的功率流并将功率流数据保存到外部文件。

由于我们感兴趣的仅是功率流本身，从而在完成功率流计算后即可提前进入下一个工况求解（Nastran允许在一次提交的分析任务中包含多个不
同载荷工况—SubCase ，它们在求解中由一个循
环控制执行）。

结构功率流分析模块的程序流程图见图
2。

根据该思路，在应用该模块分析结构功率流时，首先应将全模型提交一次频率响应分析，然后将模型分组各提交一次功率流计算分析，才能得到经过所有节点的功率流数据。

5实例分析
为了检验程序编制的正确性，使用具有解析解的简单结构进行数值计算比较。

考虑如图3所示的方形截面的悬臂梁结构，梁长L =2m ，截面宽b =10mm ，截面高h =20mm ，材料常数E =2.1×
1011Pa ，ρ=7800kg/m 3，ν=0.29，η=0.01。

在梁的自由端作用圆频率为ω的横向简谐力F =1N ，则可得
到欧拉梁模型的稳态运动方程及边界条件为
d 4
y d x 4-k 4y =0(4)
y ′′0=y ′′′0=y ″′′L =0，Q＝1+j ′′
ηEJy ′′
3=-1(5)
式中，k =ω2
ρA ′′
4
姨
，J =1bh 3
，ω是激励的圆频率。

从前式方程中可以解得梁在不同频率下的
响应复幅值y 軃′′x 。

在该问题中只存在两种内力，
即剪切力和弯矩，它们分别可以由梁的挠度曲线获得：
M ′′x =1+j ′′ηEJ y 軇″′′x ，Q ′′x =M ′′′x =1+j ′′
ηEJ y 軇′′
3′′x
(6)从而可以计算得由此两种内力引起的梁内轴向平均功率流分别为：
P 1=12
Re -Q ′′x j ωy 軇′′x ′′*
′′(7)P 2=12
Re M ′′x j ωy 軇′′′x ′′*
′′(8)
使用解析方法计算与有限元方法计算所得的前六阶固有频率比较见表1。

表1悬臂梁的固有频率单位[Hz]
Tab.1Natural frequencies of cantilever beam,unit[Hz]
图2功率流分析模块流程图，a.求解位移响
应，b.计算功率流
Fig.2Diagram on power flow analysis module,
a.solve displacement response,
b.calculate power flow
图3悬臂梁受单位简谐力激励
Fig.3Cantilever beam stimulated by harmonic
force
1
23456解析方法有限元法
4.19054.1907
26.26326.252
73.54573.463
144.12143.83
238.22237.50
355.88354.30
用两种方法分别计算梁在4.2Hz ，26.3Hz 两频率处的响应及功率流情况。

响应计算结果比较见图4，图中曲线分别表示两种方法计算得到的梁在4.2Hz 频率下的响应复数幅值的实、虚部，从中可以看出两种方法的计算结果吻合得很好。

然后，分别计算通过各个横截面的功率流—以功率从梁自由端流
第1期谢基榕等：基于有限元的功率流分析方法 (147)
向固定端为正，结果比较见图5-6，图中曲线分别给出了由剪切作用、弯曲作用引起的功率流动及其总和。

流经梁中点输向梁固定端的平均功率流计算结果比较见图7。

图4
梁在4.2Hz 处的响应图5梁在4.2Hz 处的功率流
Fig.4Beam ’s response at 4.2Hz Fig.5Power flow in beam at 4.2Hz
图6梁在26.3Hz 处的功率流图7流过梁中点的功率流
Fig.6Power flow in beam at 26.3Hz Fig.7Power flow through middle point in beam
从图7中可以看出，随着频率升高，在固有频率附近的计算结果差异增加，如26.3Hz 处的结果差异比4.2Hz 处要大，其他频率处都吻合得相当好。

图8“L ”形板在123.6Hz 下的功率流图图9“L ”形板在263.5Hz 下的功率流图
Fig.8Power flow in L-shaped plate at 123.6Hz Fig.9Power flow in L-shaped plate at 263.5Hz
148船舶力学第13卷第1期
考虑两块连接成“L ”形的平板结构，两块平板的尺寸均为1.0m ×0.5m ，将它们在长边处连接呈90度状态，将其他各边简支，平均厚度取
6.35mm ，材料常数为E =72GPa ，ρ=2710kg/m 3，η=0.01，ν=0.3。

分析此结构在集中力（距连接边及侧边均为0.2m ）激励下的功率流情况，并利用Msc/Patran 提供的结果后处理功能可以获得平板结构的功率流图及功率损耗云图，见图8-10，从这些
图中可以方便地看出板结构在振动时的功率流向、大小以及功率吸收分布等，从而为结构减振降噪设计提供重要的参考依据。

6结
论
有限元方法可以对复杂结构进行动力学分析，本文提出了在Msc/Nastran 的基础上进行二
次开发结构功率流分析模块的方法，从而克服了对复杂结构进行功率流分析的困难。

在DMAP 编程语言的基础上，将功率流分析方法以程序模块的方式予以实现并加入到Msc/Nastran 的标准求解器中，拓宽了Nastran 分析软件的应用功能。

采用解析法和本文开发的程序模块对悬臂梁结构进行的功率流分析结果表明，本文方法能够对复杂结构的功率流进行定量分析。

此外，本文还在Msc/Patran 的基础上，用PCL 语言开发了与功率流分析功能配合使用的前后处理模块，从而可以方便地生成功率流分析的数据文件和图形化显示。

参考文献：
[1]伍先俊，朱石坚.基于有限元的功率流计算及隔振系统优化设计技术研究[J].船舶力学,2005,9(4):138-145.[2]Mike Reymond.DMAP Programmer ’s Guide[K].MSC.Software Corporation,2001.
图10“L ”形板的功率损耗云图
Fig.10Power loss in L-shaped plate at
123.6Hz
第1期谢基榕等：基于有限元的功率流分析方法 (149)。