有限元理论基础及Abaqus内部实现方式研究系列42： 声学分析(1)-有限元
（原创，欢迎转载，转载请说明出处）
1 概述
本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过
（1） 基础理论
（2） 商软操作
（3） 自编程序
三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。
有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。
一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元求解器，通过自研求解器和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。
通用结构有限元CAE软件iSolver操作视频：
http://www.jishulink.com/college/video/c12884
==第42篇：声学分析==
Abaqus虽然是结构计算软件，但如今随着系统的复杂，结构与其它物理场的耦合也越来越频繁，因此，Abaqus也将自己的计算能力往其它物理场扩展，其中，热学、流体、声学和电磁学等都是Abaqus多物理场的扩展的典型应用。热学主要计算热膨胀时结构应力的变化和热传导等现象，本身和结构关系密切，所以Nastran、Abaqus等结构软件都会包括热分析功能；流体由于达索已经收购了业界两款最好的基于最近大热的LBM粒子算法的商业软件PowerFlow和Xflow，可能处于产品线的发展考虑，把基于有限元方法的Abaqus/CFD模块直接停止开发了，但保留了Abaqus中的CEL模块，专门处理流固耦合问题；而声学模块主要处理声腔内的模态分析和声传播问题；电磁模块至今不涉及，在此不讨论了。这几个多物理场的扩展每一项都需要单独的其它物理场的专业知识和结构专业的融合，也许大家都和我们一样只熟悉结构的算法，对陌生的其它专业总觉得非常难，但其实静下心来系统学习一下，其它专业的内容并不是想象的那么困难，结构有限元学好了，其它专业都是一通百通的。譬如下面的声学有限元理论公式我们就可以完全按照结构有限元的流程推导。我们在实际研发CAE软件过程中发现难的不是这些多物理场公式等的编写，难的是如何在现有的结构CAE软件基础上在前后处理和求解器层面分别加入其它分析能力，同时又保持原有框架的可扩展性和易维护性，当然，抛开多物理场，就算仅研发结构CAE软件如何在维持现有框架结构基础上在软件层面添加新的功能也是实际中我们需要花大力气研究的一个关键点，这也是我们理解的在研发上理论上的有限元和工程上的有限元分析的最大区别。
现阶段我们只讨论声学物理场，由于声学涉及的问题较多，我们仅限于和实际工作中用到的和结构强相关的部分内容，水平有限，同时对声学理论等理解的也不够透彻，可能有很多问题，大家有疑问也欢迎多讨论。我们分几篇文章来说明，本篇是声学分析（1）-有限元，讲述声学的一些基本分析算法，并着重介绍声学有限元的计算方法和融入到结构流程的关键问题，在iSolver中根据这些算法实现了声模态和稳态等声学有限元分析流程，可发现iSolver的计算结果和Abaqus没有误差，以此验证公式的正确性。
1.1 声学基本分析算法
每种学科的数值计算都有很多种，声学也不例外，在Virtual.Lab、VA One、Actran等声学商业CAE软件中，声学有限元、边界元和统计能量法是主流。
1.1.1 声学有限元
大家做过结构有限元中的振动波的问题就知道，想要模拟波的传播，那么有限元的网格就需要满足一定的要求，一般网格大小要<波长的1/20，也就是一个周期内有二十个点才能近似模拟出一条正弦曲线，下面图就是20个点组成的sin曲线，要是点少了就很难看出是正弦曲线了。
理想来说，只要计算机足够强，都可以将网格划分到波长1/20内（也许基于这个考虑，Abaqus只采用了有限元来计算声学），但实际情况是声波频率和你的分析尺度之间是很难调和的，譬如你关心的频率是300Hz，在空气中大概的波长是1m，那么1/20就是5cm，而如果你想分析一个5x5平方米房间的声传播现象，一个方向需要的网格数是5m/5cm=100，三维就是100万网格了，已经是个比较大的网格数了，所以处于实际考虑，一般声学有限元只用来处理低频问题（每个行业对低频的定义不同，我们一般认为是200Hz以下），譬如声腔内的低频模态共振情况，同时由于声音的来源很多是结构的振动和流体的脉动压力激发，对Abaqus以结构有限元为主的软件，自然而然会计算声振耦合分析问题，此时声学单元就附在结构单元上，既可以处理内音场问题，也可以处理外音场问题，当处理外音场问题时，由于外音场是无限大的，所以一般外层的声学有限单元只建1/3波长，然后外面建无反射边界或者无限元模拟无限大的情况。
1.1.2 声学边界元
即使对于低频问题，计算域大时网格数依然很多，譬如要求潜艇的几公里外的噪声分析，对这种无限大外音场问题，如果场域内部的材料不发生变化，那么边界元可以快速求解，边界元把场域看做理想的均匀材质，只需要知道源表面的振动速度，就能利用Green公式计算出任意远场的声场情况，也就是只要源靠近场域的面网格就够了，与需要把全部求解域都进行离散的有限元法相比，只离散求解域边界的边界元法更适合处理开放空间内的物理问题，它把三维问题变为二维问题求解，降低了计算复杂度。
Abaqus官方没有边界元的方法，但Abaqus集成了一个第三方的从无限元的声压计算无限域外声场的计算方法，就在帮助文档的Scripting User’s manual的9.10.11 Using infinite elements to compute and view the results of an acoustic far-field analysis。计算方法猜测还是和边界元中利用Green函数的方式类似，但我们没看懂原理，如果哪位大神了解它的理论公式，希望不吝赐教。
1.1.3 统计能量法
由于前面的基于网格的有限元和边界元只能处理低频问题，因此，在高频段完全抛弃了基于网格的计算方法，而从能量的角度来分析复杂结构在外载荷作用下的响应。它从某种程度上忽略了结构的具体细节, 同时也很好地解决了声场与结构间的耦合问题。看过别人用统计能量法来建模的过程，发现建模的随意性相对较大，同一个分析对象，统计能量法建模的实际形状的差异性很大，实在不理解为何对计算的影响可以忽略，也许这就是统计的含义，对这个方法，没深入研究过，据说Abaqus新版本里已经包括了统计能量法的模块，也没见过，如果谁用过也望告知一下入口和操作流程。
1.1.4 其它
除了低频和高频，还有中频问题，一般采用统计能量法和声学有限元结合计算。还有基于射线或者其它的不同方法，完全陌生，^.^，也只能跳过了。
1.2 声学的两个基本概念
在开始学习声学理论的时候，一直对两个概念很模糊，就是声压和声音的传播问题，在此，仅写下我们对这两个概念外行的理解。
1.2.1 声场求解物理量-声压
每个物理场都有特定的求解物理量，结构的求解物理量是各个节点的位移自由度，如果是梁壳，则加上转动自由度。声场的求解物理量是声压。声场本质上是存在于流体域中，流体在没有声场的情况下本身就有压强，譬如大家熟知的水下h米的压强p0是密ρ*g*h，奋斗者号下潜到1万米的承受的压力相当于在我们的指甲盖上站一头大象。
如果有仪器能直接测试压强，当对面有个噪声源发声时，可以发现测得的压强变为了p0+p，增加的p就是声压，声学量均为一级微量，即声压值数量级远低于介质中的静态压强。
1.2.2 声音的传播
声压是与时间和空间相关，整个区域的声压随时间变化的过程就是声音在三维空间的传播过程，声音的传播靠的是互相临近的质点在平衡位移的微笑扰动造成质点之间的挤压，质点本身并不会传播，只是一种运动形态的传播，而且介质中传播的都是小振幅声波。声音是可压缩的，否则传递声音就是无穷快，想象一个不可压缩的刚性梁左侧的力会瞬间传递到右侧。
1.3 声学有限元的理论推导
我们首先关注声学有限元，因此下面着重介绍一下声学有限元的计算方法，只不过我们换个角度，从结构有限元的角度去理解它。
1.3.1 强方程
有限元都是对偏微分方程的求解，在有限元中偏微分方程也称为强方程。
（1）对于结构某点受单位力，它的偏微分方程：
表示单位体积内的外力f是内部应力随位移的梯度，本质上还是牛顿的内外结构力平衡。
（2）对声学，Abaqus或者专门的声学软件VA One、VirtualLab等都是求解的稳态问题，即在频域求解。简谐激励下的稳态声场，声压可表述为随时间简谐振荡的量，
声学振动作为一个宏观物理现象，必须满足三个基本物理定律，即：牛顿第二定律、质量守恒定律和用于描述压强、温度与体积之间状态关系的物态方程，最终推导得到声场的偏微分方程为如下的Helmholtz方程：
上式中，
为波数，
为声波振动的圆频率，P为相应点声压幅值，Q为声源，可以表示点声源、偶极子声源或者四级子声源。具体推导可以直接看声学理论书即可，其实和流体、结构的偏微分方程推导很像。
2.3.2 弱方程
上述方程对每个场点都成立，由于实际场点是不可数的，无法采用计算机模拟，需要转化只有有限未知量的数值方程，第一步就是转换为弱方程。
（1）对结构来说：
Step1：转换为弱方程
在每个网格内既然偏微分方程成立，那么乘以任意一个矢量函数w对网格体积积分依然成立：
上式即是采用伽辽金方法获得偏微分方程的加权余量形式，w称为权重函数。注意，原有偏微分乘以权重函数要使得结果是标量，因为原有偏微分方程是一个向量，所以这边权重也是一个同等维数的向量，上述为点乘。
Step2：降阶
等参形式下，任意点的位移u可以表示网格节点位移ui的线性组合。
上述实际是三个方向的三个方程，分解出来应该是下面的式子：
应变是u的一次导数，因为ui与坐标无关，那么带入应变表达式后可以转换到一次偏导就放到形函数上：
上式不严谨，只表达应变中与位移自变量相关的是形函数的偏导，如果三维情况真实就是：
u的一次导数是将对xyz偏导放到了形函数上，但对二次导数就无法这么做了，以线性材料为例，应力和应变是线性关系，
就是一个对xyz的二次导数，如果再次把对位移的偏导转到中，此时变为
如果是一个线性一次单元，那么上式恒等于0，和施加的外力f就不平衡了。
所以u的二次导数需要降维，将二次导数转化为一次导数。利用的就是大家熟知的gauss定理，该物理量的散度在体积V内的体积分=物理量在包围面S的面积分。散度即偏导，高斯定理的意义就是降了一次偏导（散度）。
其中n从系统域指向外面。
则方程（1）转换为
为了更好的说明上面的含义，w是任意的函数，在这儿，不妨取w为u，且内力和外力项分别放等式左右两边，得到
即为应变的变分，则左侧相当于位移变化
，
即为每个点的应变能的变化，
为整个系统的应变能变化，而右侧即为外力做的虚功。所以上式即为虚功原理，在物理上可解释能量守恒原理，即在某一个时刻点，假定在外力作用下有个虚拟的位移，那么外力在虚拟位移下做的虚功=内部应变能的变化相同。
（2）对声学：
Step1：转换为弱方程
在每个网格内既然Helmoheltz偏微分方程成立，那么乘以任意一个标量函数w对网格体积积分依然成立：
Step2：降阶
等参形式下，任意点的声压P可以表示网格节点声压Pi的线性组合。
声学的形函数可以和结构取成完全一致的，最常用的就是线性形函数，此时
就恒为0。所以也需要用高斯定理降阶。
等式（2）就变为
上式只要将
和
看成(1)式中的
和f就可得到
1.3.3 有限元近似
强方程转化为弱方程后，转化只有有限未知量的数值方程，第二步就是有限元近似，我们选择一些权重函数w，这些权重函数不是在系统域内都是任意的，而只是有限个点任意，但满足两点：
严格的系统域内积分可以表达为另一种形式的函数，这些函数的未知量只有这些有限点上的物理量，这样可以采用计算机求解。
当这些有限点的数目足够密时，严格的系统域内积分就可以近似为另一种形式的函数了。
有限元就是其中的一种取任意点的表达形式。
按照有限元的近似，可以将系统域划分为有限的单元，可以认为原有系统域就和有限的单元的体积累加起来近似。
同时，也认为严格的弱方程对原有系统域的积分函数近似为有限的单元的体积域积分函数累加：
（1）对结构得到
其中，上面的边界不是每个单元都有值，只有边界在系统域边界S上的才有值。
权重函数w在有限元理论中可由为节点插值得到，且插值函数和u相同
可得：
代入上式变为：
其中
只是简化形式，实际应该是上面应变和位移的关系矩阵[B]。
因为
为任意函数，所以必须满足：
此时中的t、
和f仅包括未知量ui，即上式为一个只包括ui的代数方程组，利用线性或者非线性方程组的一般求法即可求出。
（2）对声学，将声学弱方程近似为有限元方程后得到
令权重函数采用同样的形函数插值
代入上式得到
由于
是任意的，所以得到
其中
和P仅包括未知量节点声压Pi，上式即是一个只包括Pi的代数方程组。
最后再聊一下声学边界情况，由声学边界条件：
同时，Abaqus和iSolver的声学原始方程整体都除了一个
，这样边界就只有加速度了。即声学边界上的载荷可以由从边界指向声场内部的法向加速度（inward volume acceleration）确定，所以在Abaqus或者iSolver中，都采用了该法向加速度作为声学分析的边界载荷。
2.4 基于iSolver结构流程的声学有限元分析实现
一般情况声学方程都是线性代数方程组，不需要迭代就可以求出，但因为iSolver求解器已有增量迭代法的结构求解流程，我们程序实现中还是按迭代来求解，这样我们只要加入了声学单元，同时求解声学单元的刚度阵、质量阵及非平衡力，只不过一次迭代就收敛结束了。
为和Abaqus的声学方程一致，在iSolver实际程序代码中，我们是将原始方程*负号，表示为：
1.5 声学有限元的卡车噪声模型验证
1.5.1 模型介绍
该模型分析卡车受到地面对轮胎的激励导致的噪声分析。该模型从郑钧Adam老师的Abaqus線性動力&噪音分析詳解(理論及實作) Abaqus 線性 動 力 & 噪音分析詳解 ( 理論及實作 ) 视频教程 _ 培训课程 - 技术邻 (jishulink.com) 下载获得。
由于卡车为对称模型，所以只用建一半模型，声音到了对称面将全反射，但和结构分析不同，声学有限元边界默认全反射，所以在此对称面无需任何额外设置。腔体内为空气材料，全域划分为四面体声学网格，供11681个单元。
1.5.2 模态分析
在iSolver中创建模态分析步，设置为0.1Hz-200Hz分析区域，计算15阶模态。
1.5.2.1 iSolver结果
iSolver计算完毕在该区域内只有三阶频率，分别为：
第一阶76.6Hz。
第二阶134.2Hz
和第三阶168.5Hz
1.5.2.2 Abaqus结果
在Abaqus做通用的设置，得到的模态个数和频率完全一致：
1.5.3 谐响应分析
在轮胎两个位置所在的节点设置加速度载荷，采用基于直接法的稳态动力学分析步，设置50Hz-200Hz内20个频率点。
1.5.3.1 iSolver结果
采用iSolver计算，计算完毕查看车位头部的声压。
在iSolver中绘制声压随频率变化的曲线结果如下：
1.5.3.2 Abaqus结果
Abaqus计算后选择同样节点得到曲线如下：
可发现在76Hz和134Hz两个频率点附近Abaqus和iSolver均有共振峰，但第三阶模态均没激发出来。从两个软件的分析结果来看，iSolver声学模态及稳态分析的结果和Abaqus精度完全一致。
1.6 视频讲解和操作验证演示
如果觉得上面的文字太复杂，也可以看一下视频的简要讲解和软件演示，地址如下：
https://www.jishulink.com/college/video/c12884 20理论系列文章42-声学分析(1)-有限元
1.7 以往的系列文章
1.7.1 ========第一阶段========
第一篇：S4壳单元刚度矩阵研究。
http://www.jishulink.com/content/post/338859
第二篇：S4壳单元质量矩阵研究。
http://www.jishulink.com/content/post/343905
第三篇：S4壳单元的剪切自锁和沙漏控制。
http://www.jishulink.com/content/post/350865
第四篇：非线性问题的求解。
http://www.jishulink.com/content/post/360565
第五篇：单元正确性验证。
https://www.jishulink.com/content/post/373743
第六篇：General梁单元的刚度矩阵。
https://www.jishulink.com/content/post/403932
第七篇：C3D8六面体单元的刚度矩阵。
https://www.jishulink.com/content/post/430177
第八篇：UMAT用户子程序开发步骤。
https://www.jishulink.com/content/post/432848
第九篇：编写线性UMAT Step By Step。
http://www.jishulink.com/content/post/440874
第十篇：耦合约束（Coupling constraints）的研究。
https://www.jishulink.com/content/post/531029
1.7.2 ========第二阶段========
第十一篇：自主CAE开发实战经验第一阶段总结。
http://www.jishulink.com/content/post/532475
第十二篇：几何梁单元的刚度矩阵。
http://www.jishulink.com/content/post/534362
第十三篇：显式和隐式的区别。
http://www.jishulink.com/content/post/537154
第十四篇：壳的应力方向。
https://www.jishulink.com/content/post/1189260
第十五篇：壳的剪切应力。
https://www.jishulink.com/content/post/1191641
第十六篇：Part、Instance与Assembly。
https://www.jishulink.com/content/post/1195061
第十七篇：几何非线性的物理含义。
https://www.jishulink.com/content/post/1198459
第十八篇：几何非线性的应变。
https://www.jishulink.com/content/post/1201375
第十九篇：Abaqus几何非线性的设置和后台。
http://www.jishulink.com/content/post/1203064
第二十篇：UEL用户子程序开发步骤。
https://www.jishulink.com/content/post/1204261
1.7.3 ========第三阶段========
第二十一篇：自主CAE开发实战经验第二阶段总结。
https://www.jishulink.com/content/post/1204970
第二十二篇：几何非线性的刚度矩阵求解。
http://www.jishulink.com/content/post/1254435
第二十三篇：编写简单面内拉伸问题UEL Step By Step
http://www.jishulink.com/content/post/1256835
第二十四篇：显式求解Step By Step。
https://www.jishulink.com/content/post/1261165
第二十五篇：显式分析的稳定时间增量。
http://www.jishulink.com/content/post/1263601
第二十六篇：编写线性VUMAT Step By Step。
https://www.jishulink.com/content/post/1266640
第二十七篇：Abaqus内部计算和显示的应变。
https://www.jishulink.com/content/post/1273788
第二十八篇：几何非线性的T.L.和U.L.描述方法
https://www.jishulink.com/content/post/1282956
第二十九篇：几何非线性的T.L.和U.L.转换关系
https://www.jishulink.com/content/post/1286065
第三十篇：谐响应分析原理
https://www.jishulink.com/content/post/1290151
1.7.4 ========第四阶段========
第三十一篇：自主CAE开发实战经验第三阶段总结
https://www.jishulink.com/content/post/1294824
第三十二篇：谐响应分析算法
https://www.jishulink.com/content/post/1299983
第三十三篇：线性瞬态动力学
https://www.jishulink.com/content/post/1302074
第三十四篇：非线性瞬态分析
https://www.jishulink.com/content/post/1787283
第三十五篇：接触求解算法
https://www.jishulink.com/content/post/1792869
第三十六篇：DLOAD用户子程序开发步骤
https://www.jishulink.com/content/post/1826803
第三十七篇：梁单元差异（1）-理论基础
https://jishulink.com/content/post/1872208
第三十八篇：梁单元差异（2）-梁截面方向
https://www.jishulink.com/content/post/1874628
第三十九篇：梁单元差异（3）-剪力和弯矩
https://www.jishulink.com/content/post/1876013
第四十篇：梁单元差异（4）-形心、剪心和偏置
https://www.jishulink.com/post/1888000
1.7.5 ========第五阶段========
第四十一篇：自主CAE开发实战经验第四阶段总结
https://www.jishulink.com/post/1905963
