声场求解器开发之(3)---C 开发BEM求解器

本文主要介绍基于BEM（边界单元方法）方法开发用于计算声学求解器的基本知识。

1. 概要

在数值仿真领域，有限元方法一直是首选，但在一些细分领，其他方法可能更有效，比如有限体积法（FVM）对流体，边界元法（Boundary Element Method）对声，电磁散射。

以声场为例，声音向外传播，作用区间是一个典型的无限范围，如果使用FEM，需要对整个区域划分网格，尤其在三维空间，随着计算范围的增大，网格数目会急剧增大。而采用BEM，只需要求解声场边界上的数值，大大降低了计算量，提高了计算效率。

既然BEM效率这么高，为什么使用范围不如FEM广泛呢？主要两点原因：

1. FEM 是一种纯数值解法，而BEM是半数值解，也就是说要以解析基本解为前提条件

2. 因为最后形成的线性方程组为非对称满秩矩阵，限制了工程上的求解规模。桌面单机几千自由度就已无法解出结果，虽然通过FMM（快速多极子方法）解决了这个问题，但因为FMM算法本身实现也有一定难度，因而没有FEM应用广泛。统计显示当模型自由度在15000以下时，传统的BEM效率更高，而大于15000时，传统BEM计算量会急剧上升，此时FMM的计算量仍然呈线性增长。

BEM 基于解析解，在处理某些特定领域（声场，电磁场，连续介质弹性力等无限域问题）具有 精度高，降低维度等特点，同时通过与FEM相结合，能够综合两者的优势，提高计算效率和精度。

2. 理论

前面介绍过声场的控制方程为赫姆霍兹方程，边界元方法的第一步是要将PDE转换为边界积分方程，对于三维空间的方程，需要使用格林函数将体积分转换为面积分。

通常流程是：

1. 利用算子的基本解作为权函数，按加权余量格式得到区域上的积分方程；

2. 利用高斯公式（格林函数）建立区域内积分和边界积分的关系，从而得到区域内任意一点的通解变量表示的积分表达式；（利用格林函数可以将三维体积分转为二维面积分）

3. 将基本解的奇异点P趋于边界点p，得到边界积分方程。

这个过程简称为用加权余量法建立边界积分方程。

BEM求解流程与FEM基本类似，划分网格，建立线性方程组，求解方程，不同的只是建立的本构方程不同。

因为BEM求解区域只在边界上，所以三维问题只需要划分面网格，二维问题划分线网格，使前处理工作大为简化。

BEM的计算流程与FEM相似：

1. 生成网格（三维使用面单元，二维计算使用线单元）

2. 计算单元影响矩阵

3. 组装总体影响矩阵

4. 求解线性方程组

5. 根据边界上的值计算任一点的声场压力

不失一般性，以线性四边形单元求解空间声场问题。

单元定义可以从FEM的四边形单元派生，因此可以直接使用FEM中等参单元定义的形函数， 再加上声场计算所需的函数与变量。针对每个单元需要计算H和G矩阵，这两个矩阵后面用来组装整体矩阵。为了求解精度，三维实体划分网格的时候，四边形边长需要控制在波长的一定比例之内。

1. 声场线性四边形单元定义如下：

单元函数说明：

CalGeMatrix() 和 CalHeMatrix()分别用来计算 G矩阵和H矩阵。关于矩阵介绍参考附录 参考1。

2. 材料属性：

材料可以直接使用FEM定义的部分

1. 角频率 = 2*pi*v (v为频率 单位 赫兹)；

2. 声在介质中的传播速度；

3. 介质的密度；

3. 整体刚度组装流程：

4. 边界条件带入总体矩阵，构建线性方程组

边界条件包括压力，速度以及声阻

5. 结果计算

求解得到边界上的声场压力和速度之后，就可以直接利用表达式计算得到空间中任意一点的压力和速度。

由于最后的HG矩阵为满秩矩阵，在本人I7四核 8G PC机上单元达到5000时，基本上已经算不出结果。

参考：

1. STRUCTURE-ACOUSTIC ANALYSIS USING BEM/FEM IMPLEMENTATION IN MATLAB. FREDRIK HOLMSTRÖM

2. The Boundary Element Method in Acoustics. Stephen Kirkup

3. BEM-FEM Acoustic-Structure Interaction For Modeling and Analysis of Spacecraft Structures Subject to Acoustic Excitation.   Harijono Djojodihardjo

声明：原创文章，欢迎留言与我讨论，如需转载留言