基于iSolver的声学有限元分析（2）| 边界元
1 概述
本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过
基础理论 商软操作 自编程序
三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。
有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。
一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元软件iSolver，通过自研CAE软件和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。
2 声学分析
Abaqus虽然是结构计算软件，但如今随着系统的复杂，结构与其它物理场的耦合也越来越频繁，因此，Abaqus也将自己的计算能力往其它物理场扩展，其中，热学、声学、流体和电磁学等都是Abaqus多物理场的扩展的典型应用。
热学主要计算热膨胀时结构应力的变化和热传导等现象，本身和结构关系密切，所以Nastran、Abaqus等结构软件都会包括热分析功能；流体由于达索已经收购了业界两款最好的基于最近大热的LBM粒子算法的商业软件PowerFlow和Xflow，可能处于产品线的发展考虑，把基于有限元方法的Abaqus/CFD模块直接停止开发了，但保留了Abaqus中的CEL模块，专门处理流固耦合问题；
而声学模块主要处理声腔内的模态分析和声传播问题；电磁模块至今不涉及，在此不讨论了。这几个多物理场的扩展每一项都需要单独的其它物理场的专业知识和结构专业的融合。
也许大家都和我们一样只熟悉结构的算法，对陌生的其它专业总觉得非常难，但其实静下心来系统学习一下，其它专业的内容并不是想象的那么困难，结构有限元学好了，其它专业都是一通百通的。
譬如下面的声学有限元理论公式我们就可以完全按照结构有限元的流程推导。我们在实际研发CAE软件过程中发现难的不是这些多物理场公式等的编写，难的是如何在现有的结构CAE软件基础上在前后处理和求解器层面分别加入其它分析能力，同时又保持原有框架的可扩展性和易维护性，当然，抛开多物理场，就算仅研发结构CAE软件如何在维持现有框架结构基础上在软件层面添加新的功能也是实际中我们需要花大力气研究的一个关键点，这也是我们理解的在研发理论上的有限元和工程上的有限元分析的最大区别。
现阶段我们只讨论声学物理场，由于声学涉及的问题较多，我们仅限于和实际工作中用到的和结构强相关的部分内容，水平有限，同时对声学理论等理解的也不够透彻，可能有很多问题，大家有疑问也欢迎多讨论。分几篇文章来说明，本篇是声学分析（2）-边界元。
由于声腔区域很多都是均匀的，这个区域很容易只要表达出边界情况就能确定解，所以边界元在声学中经常用到，iSolver仅仅聚焦有限元，对于声学有限元，我们团队自己将结构有限元改造，加入了基于有限元的声学分析模块，但对于声学边界元，我们没有做实际的开发，也不打算在边界元领域深耕，只是希望和其它团队合作在iSolver上二次开发。
当然，很多人连自主CAE软件直接操作都懒得试一下，更不要说愿意投入大把时间基于自主CAE软件二次开发，但只要有人有意愿在iSolver平台上二次开发，我们都会鼎力支持，现在的iSolver平台已完全支持基于python和C++的二次开发。
另一个团队在iSolver前后处理平台上二次开发完成声学边界元模块，因为涉及界面的边界元整体流程和底层的数据接口、声振耦合传递等，所以我们对声学边界元略有了解，但比起真正做边界元的还是差很多，仅记录下我们理解的边界元公式和开发难点。
2.1 声学边界元方程推导
声学边界元的标准推导可以直接看理论的书，我们只懂有限元，所以照着边界元的理论书尽量按上一篇有限元的强方程-->弱方程-->有限元近似的流程重新简要推导了一下。
2.1.1 强方程声学边界元方程的推导最主要的是如何把一个体内域V的计算直接转为只要知道边界上的物理量就行。还是从Helmholtz声学波动方程出发：
这个方程有个特解，就是当一个源强为1的点源位于r0点时：
此时显然这个P围绕r0点是球对称分布，也就是只与|r-r_0 |有关，将上式用球坐标表示后，可以求出这个特解为Green函数
即上式的G满足：
2.1.2 积分方程
一般，边界元解决的工程问题都是结构表面振动或者流体表面脉动压力导致声音在外空间均匀介质传播的问题，譬如汽车盖在空气中由于发动机振动导致的噪声、潜艇在无限水域的辐射噪声。此时声源只在表面，辐射区域没有声源，对V内的任意一点。
其中 为辐射体域内任意一个场点位置。
再假设这个辐射面是封闭的，这类问题也称为直接边界元问题。设这个辐射域由两个面组成，一个是结构表面Sa，另一个应该是无穷大的球面SR2。
因为声源只在表面，所以（2）式变为：
其中表示在r处有源强为1的声源， 还是为辐射体域内任意一个场点位置。用一个小球SR1包围r点的声源。在Sa（去掉SR1相交的部分，当SR1无穷小时可忽略）、SR1（去掉和Sa相交的部分）和SR2包围的辐射域Va内。
得到
对上述（1）（2）在体积Va内积分得到：
由上一章高斯定理：
体积分可以转换为面积分，得到：
约掉∇P*∇G项
当SR1足够小时，在光滑表面SR1小球和Sa相交成半球，所以右边
如果r位于角上，则需要修正，否则结果和商软将差异较大。最终转换为只有边界上的积分：
粗看这个方程神奇的把k平方去掉了，但实际只不过转换到了G函数中。
2.1.3 边界元近似
上述方程个人觉得不能称为有限元中的弱方程，因为不是像有限元一样强方程乘以权重函数来形成积分形式的，但这个积分形式的离散可以和有限元类似，最终目点还是但满足以下两点：
严格的系统域内积分可以表达为另一种形式的函数，这些函数的未知量只有这些有限点上的物理量，这样可以采用计算机求解。 当这些有限点的数目足够密时，严格的系统域内积分就可以近似为另一种形式的函数了。
将系统边界Sa划分为有限的边界单元，可以认为原有边界域就和有限的边界单元的面积累加起来近似。
同时，也认为严格的积分方程对原有系统域的积分函数近似为有限的单元的体积域积分函数累加：
得到：
我们一般处理频域问题，所以可以设
由声学边界条件：
Vn是法向速度，从声腔指向振动源，方向如上图所示。
在实际编程中，如果不是做纯有限元，而是结构或者流体导致的外场问题，那么法向尽量取成从结构或者流体指向声腔，因为一般结构或者流体体单元时，体上的面默认都是指向体的外面的。
我们还是以符号n为指向流体域的方向，得到：
即：
结构或者流体导致的外场问题，一般都已知法向速度Vn，求声压P（r）。划分为网格后，设网格点上的坐标为rn，只要将网格点rn对应的声压求出来就行。
上述还需解决单元内未知量P(ra)，和有限元类似，对一个单元内的声压P(ra)采用边界元节点表示，等参形式下，任意点的声压P可以表示网格节点声压Pi的线性组合。
这样，单元内的未知量就直接转换为了节点的声压未知量组合了。
将r=rb，b=1…N（N为边界上节点总数），遍历，得到N个线性方程组，且包括N个P(rb)未知量。
上式直接求解可得Pi。
2.2 基于iSolver结构流程的声学边界元的程序实现的关键
边界元的原理相对有限元来说并不复杂，但要做出一个能带操作界面且精度效率可靠的工程软件还是非常难的。合作团队和平台都投入了大量时间来共同实现基于iSolver的声学边界元的代码开发，我们理解边界元分析存在以下几个关键技术：
2.2.1 大模型的存储和计算问题
边界元第一个问题是求解精度和内存消耗，其实抛开内存谈效率和精度没有意义，由上述最终的方程可知，待求的P(r)和所有单元的（ ）有关，导致得到的P(r)前系数阵是个满阵，满阵的存储内存和计算都是个问题。
在iSolver中实测1G内存大概只能存下1500个边界单元，40G内存大概存下1w个边界单元，基本上内存开销和单元个数成平方比，对于计算效率，虽然边界元涉及的矩阵均为非稀疏矩阵，但实际上我们发现现在成熟矩阵数学库一般已经能兼容稀疏阵和非稀疏阵，如果工程使用的话直接调用成熟库就行。
当然，有能力也可自己编写更快更适合自己实际情况的矩阵求解。现在商用声学软件都采用FMM快速多级子通过减少节点间的相互关系来减少内存和计算加速，缺点是也会损失些精度，而且算法复杂编程实现和维护相对也要花费更多时间。
2.2.2 声学边界元和结构有限元耦合分析问题
第二个问题是前后处理的流程：工程实际的噪声源一般都是结构CAE或者流体CFD计算得到的边界振速，对于声学量会影响噪声源的情况，需要双向耦合，而对于影响较小的，可以直接将这些CAE软件得到的边界传递到边界元的网格上，这些都涉及到两套网格之间的匹配问题，网格不同需要通过形函数等插值保持力和力矩平衡。
实际工程问题可以不单独构建边界元网格，而是直接在结构或者流体的CAE网格抽取耦合面的面网格作为边界元网格，在Virtual.Lab等声学边界元商软中很多工程就是这么做的，缺点是CAE网格一般比较密，导致抽取的边界元网格相对较多，容易导致前面大模型导致的存储和计算问题。
就算是通过抽取表面的简单方式来耦合，iSolver在实际前后处理编程中也为如何操作耦合考虑了很久，毕竟Abaqus等结构商软没有边界元的现成流程可参考，同时Virtual.Lab也是直接导入CAE后处理开始的，最终iSolver前后处理的流程实现边界元有3个关键步骤：
直接从一个part模型上抽取表面形成边界元 
类似稳态动力学的声学边界元的Step创建 
采用子模型的设置中来关联边界元和结构后处理数据。 
2.2.3 融入到有限元分析流程中的边界元求解过程
上面的是前后处理的关键流程，还有求解器的问题：对边界元求解流程如何融入到有限元分析流程也是一个难点。上一章讲过声学有限元只要加入声学单元，求出类似的刚度阵和非平衡力，就很容易嵌入到基于增量迭代法的有限元结构流程中，但边界元实际上融入这个流程还有相当的困难.
按照最终的方程来说，我们可以把P(r)前的系数阵当成刚度阵，然后也可以采用迭代法来求非平衡力，正常来说也是一次平衡，但边界元基本不这么做，我们理解困难点在于全局刚度阵的组装，有限元中由于节点只与跟它相连的单元节点影响，可以先求出单元刚度阵得到该单元内节点之间关系，然后只需要对结构单元一次遍历就能把所有单元间节点关系累加得到全局刚度阵，也就是节点和所有其它节点之间的关系，但对边界元来说，因为每个单元节点未知量都和其它所有单元相关，所以事实上没有单元刚度阵一说，需要对单元两次遍历才能累加所有的节点间关系，也就无法用有限元稀疏阵的组装方式。
这看起来是个很小的改动，但当一个大型程序的流程固定后，看似一点小改动在代码上也很难共存，当然，相信Nastran和Abaqus等商软绝对有能力来改造，只是最终没实现而已，不清楚具体原因，也许是应用领域比较狭窄和编程成本的平衡考虑吧，一种简单方式是类似Nastran的Glue或者Abaqus的Tie一样单独开辟一个节点间关系的流程。
边界元合作团队在实现边界元过程中采用的是另一个单独的exe来完成边界元求解，求解输入输出数据文件按iSolver的以往的规范，这样可以同步采用iSolver的前后处理，当然，就算融入到以结构为主的前后处理软件中，边界元单元、属性、材料等的设置能否和结构一体化也是大家需要考虑的问题。不管如何，让用户感觉一体化是最重要的，类似Abaqus的standard和explicit求解器融入到Abaqus/CAE的前后处理一样。
2.3 直接边界元脉动球外场辐射模型校核
无论什么计算模块集成到iSolver中，我们要求和自己的求解器iSolver一样，都首先必须保证精度。我们仅以一个简单算例考核一下集成到iSolver的声学边界元求解器精度。
2.3.1 模型介绍
以一个半径为0.5m的脉动球为例来分别验证边界元计算结果。脉动球网格划分网格如下所示，网格数为900，网格类型均为四边形网格；声学介质为水，密度为1000 Kg/m3，体积模量为2.18e+09 Pa。
2.3.2 边界元分析步设置
iSolver中设置为边界元分析步，且外场计算，频率点设置为100-900Hz，取9个频率点。
创建外场中圆柱场点，圆柱壳半径5.0，长度10.0。
设置球表面法向振速为2m/s。
2.3.3 自主边界元结果
采用边界元求解后导入iSolver后处理查看云图，选择Magnitude幅值，iSolver中得到900Hz的声压幅值云图如下，在脉动球表面声压幅值=2.92e6Pa：
2.3.4 和VirtualLab对比
将iSolver模型导出为bdf，然后导入声学商业软件Virtual.Lab.Acoustics中，并施加同样的振动速度，Virtual.Lab中脉动球900Hz表面声压结果为2.93e6和2.97e6之间。
2.3.5 和理论对比
脉动球外场声压解析表达式为：
上式中， 为声学介质密度， 为声速，k=ω/c， a为脉动球半径， r为场点到球心的距离。当r=a时，可得在脉动球表面在900Hz声压理论幅值= 2.92e6；
2.4 总结
本章我们首先推到了声学边界元的数值方程，然后合作团队通过iSolver上的二次开发将边界元流程融入到了以结构有限元为主的iSolver前后处理和求解器中，通过典型算例可得边界元数值结果和边界元商软Virtual.Lab及理论基本一致。
其实在iSolver中还实现了结构振动和声学边界元的耦合的操作流程，并且结构有限元结果和边界元载荷的后台数据实现了自动流转，这才是真正实际工程的应用场景。最后，感谢我们的合作团队的辛苦努力，自主CAE软件平台不好用，二次开发更是问题多多，理解我们的有限元的数据结构和接口是绝对需要花很多时间的，当然，数据结构都是我们自己编的，所以可以比商软解释的更清楚，而且起码没有原则上改不了的地方，这也是我们自主CAE的优势之一。
