CN104239676A - 在宽频范围对系统声发射与传播建模的方法和计算机产品 - Google Patents

Abstract

本发明涉及在宽频率范围内预测声源发射及其在周围介质的传播的处理，包括源和周围介质的系统由小且简单的元素e表示。对于每个元素e和每个频率f_i，先验误差估计器使描述近似在频率f_i的声场的多项式次数的参数P_e,i与元素关联，确定对应于估计器计算的在频率范围内的最大P_e,i参数的元素的参数P_e,max。对于每个元素e，使用所述参数P_e,max确定基本矩阵K_e,max和M_e,max，它们描述元素分别对系统刚度和质量的贡献。对于每个频率f_i且对于每个元素e，参数P_e,i用来确定近似声场的多项式次数，对于所有元素从矩阵K_e,max和M_e,max提取基本矩阵K_e,i和M_e,i，其组装成分别表示系统的刚度和质量的全局矩阵K_i和M_i，基于全局矩阵K_i和M_i建立全局矩阵系统Z_i，及最终使用线性求解器对全局矩阵系统求解。

Description

在宽频范围对系统声发射与传播建模的方法和计算机产品

技术领域

本发明涉及预测周围介质内声源的发射及其传播的领域，源和介质形成系统。应当在宽频率范围内执行预测。系统可以是任何一种系统，基本上为机械系统，如车辆或电子设备。

背景技术

噪声对健康的影响越来越多地被认识到，并且成为公共健康问题。的确，暴露于升高的声级公认是极为危险的，并且可能造成听力损伤、高血压和睡眠障碍等。最显著的风险是由车辆和飞机噪声、长期暴露于大声的音乐以及工业噪声引起的。

这些考虑已经致使噪声降低成为对于当今的制造商的主流问题。例如，在汽车行业，声音发射是当今汽车设计中成熟的规范。首先，客户为了他们自己的舒适希望获得更加安静的产品。其次，授权噪声限值对于汽车、卡车和公共汽车非常严格，此外还有其他源。

目前，噪声降低仍然受到缺乏高效的预测建模工具来模拟机械系统的声音发射的限制。改进计算声模拟对于优化产品设计是必不可少的，同时避免后期的和昂贵的物理测试。图1呈现了在给定的频率下汽车腔体内部所计算的声压，其中区域越暗，则腔体内声压越高。

根据定义，声音是作为由听觉范围内的频率组成的、通过固体、液体或气体传送的压力的振荡的机械波。人类可感知的声音具有从大约20Hz至20,000Hz的频率范围。在空气中，在标准温度和压力下，声波的相应的波长范围为从17m至17mm。如果在系统的所有点处的声压和声速已知，则在给定系统中的声场适当地被定义。

在介质中的声波传播是复杂的现象，为了描述该现象，科学家基于若干近似已经构建出物理模型。由一组偏微分方程来近似波传播的模拟。这些方程的解析解没有简单解，并且必须使用计算方法来近似声场的解。

数值法的所有主要的类别已经应用于声学，包括有限元方法(FEM)，有限差分法(FDM)，间断伽辽金方法(DGM)和边界元法(BEM)。可以基于时域或频域求解器的使用对这些方法进行分类。频域求解器更适合永久物理条件(引擎运行)，而时域求解器需要研究瞬态应用(关门)。此外，可以区分需要完全呈现传播的流体区域的域方法以及依赖于初始问题的边界的积分表示的边界方法。有限元方法和间断伽辽金方法自然地处理将在下文定义的非结构网格，并且更适合于复杂的现实生活中的工程问题。FEM通常适合于频域的应用而DGM适合于瞬态应用，但是DGM也可以在频域背景中使用。

FEM仍然是工业部门中最流行的方法。其是用于求解偏微分方程的数值方法。标准的FEM过程的特征可以由以下三个主要步骤来描述：

1)将给定的机械系统划分为许多小的、非重叠的和简单的“元素”。这些简单元素或子域的集和被称为“网格”。网格描述几何特征，并且随后会包括定义波发射和传播的模拟的问题的性质和边界条件。创建网格的处理被称为网格化(meshing)。商业FEM和DGM解决方案通常需要和处置不同形状的元素(三角形，四边形，六面体，四面体，棱柱，棱锥...)。图2和图3作为图示给出了同一汽车腔体的两种网格。

2)然后针对这些元素中的每个元素来构建分段多项式近似(基于形函数(shape function))以近似感兴趣的场(在我们的情况下为压力或声速)。图4公开了形函数的16种示例。形函数对场的贡献是近似问题的未知数。这些离散的未知数通过牵涉耗时的求面积规则的一组基本方程在每个子域上链接。所有这些局部贡献(基本矩阵)然后重组成方程组的全局系统(系统矩阵)。所有这些处理被称为系统的组装。

在图4中，这些形函数中的一些相当简单，看起来像一个角升高的平面正方形。更难以描述地，其他包括凸起和凹陷，表示分布在元素域中的局部最小值和最大值。本领域的技术人员熟悉这些形函数。这些形函数建立了将在下文中解释的分层基础。

方程的求解体现在牵涉所有形函数对全局问题的贡献的一组线性方程的求解。在该步骤使用稀疏线性求解器。在每个频率组装并且求解的系统矩阵表示为Z(f)。该矩阵的维度直接与系统中形函数的总数相关，系统中形函数的总数通过网格中元素的数目乘以每个元素内的多项式形函数的数目来给出。在大部分情况下，矩阵采取以下形式：

Z(f)＝K-(2πf)²M+C(f)

其中，K是所谓的刚度矩阵，M是所谓的质量矩阵，以及C(f)包括通常来自边界条件的所有其他项。重要的是注意，质量矩阵和刚度矩阵通常与频率无关，因此它们只需要计算一次并且不需要在每个频率计算。这是在标准的FEM商业解决方案中通常所做的。

在网格化期间定义的每个简单元素能够通过以下两个主要特征来描述：

-元素的维度h，其是几何参数，

-用于对解进行近似的多项式的多项式次数P(也称为给定元素的内插的阶)。

这些参数h和P的选择是在声学仿真中要考虑的重要方面。事实上，当在常需要高阶(高P)方法和多频率解的频域中进行求解时，重要的是提出高效的策略以避免在每个频率重复地组装方程的代数系统。在这些情况下，该组装的步骤的确会是非常计算密集的。

可以通过与所研究的声音的频率的比较来完成h的选择。当在频域中进行求解时，必须独立地对感兴趣的每个单独的频率进行求解。由用户以升序存储的列表的形式来提供频率的列表。频率的总数N_f取决于应用，并且其可以是大的，因为用户可能希望以精细频率增量覆盖完整的音频频率范围。在标准低阶FEM中，基于频率来设计网格化操作以便保证以感兴趣的频率传播的声平面波具有足够的分辨率。该“经验法则”是近似，因为解不是先验已知的，并且会比仅平面波的情况更复杂。对于线性有限元，可以由以下方程来近似h：

$h=\frac{c_0}{6f}$

其中c₀是声音在传播介质中的速度，f是待研究的频率，以及h是网格化操作所需的典型网格维度。

该式表示低频率(大h)的要求低于高频率。在图2和图3中示出了这一点，其中相应地，低频率需要稀疏网格(图2)，而高频率需要精细网格(图3)。显然，与稀疏网格相关联的数值成本比与细网格相关联的数值成本小几个数量级，因为元素和待求解的相应的方程的数目会小得多。然而出于实际原因，主要是因为网格化操作是繁琐的并且不能自动化，用户常使用单个网格(精细网格)来计算整个频率范围这意味着在性能方面的巨大损失，因为低频率通常会被过度求解。

可以基于局部被呈现的场的复杂度的一些知识来局部地适应P的选择。例如，如果场是平滑的并且有很多振荡，则将局部地需要高阶近似(即，牵涉许多形函数)。找到正确数目的形函数以应用于每个单独的元素是至关紧要的，因为方法的计算成本(组装和求解过程)与牵涉的形状的数目成正比。

存在两种固定元素内的阶的方法：通过先验误差估计器或通过后验误差估计器。如由其名称所表示，估计器估计给定解的精确度，并且帮助“建模工程师”确定他们的参数P的集合。后验估计器在方程求解后细化参数P，而先验求解器在方程求解之前确定参数P。

先验误差估计器：

在实际感兴趣的应用中，在计算之前不能够猜测确切的解。然而，复杂解基于在计算前已知的控制方程的简单基本解(如平面波解)的线性组合。可以得出“经验法则”来固定阶，正如标准FEM固定网格尺寸。例如在低频率，声波波长通常较大，因此能够使用低多项式阶(较少形函数)。在高频率，波通常发生更多振荡，并且应当使用较高阶。

后验误差估计器：

在该情况下，通常以低阶计算第一廉价解。基于该第一粗略解的后验分析(查看例如初始运算符的残差或解的导数的行为)，每个元素内的阶通常会升高。计算另一解，以此类推，直到后验误差处理判断解满足所有位置。这被称为P自适应法，其中生成具有递增复杂度的收敛解的序列以对单一频率进行求解。P自适应法还没有被广泛地这样使用，它们能被看做是更一般的hP自适应法的特定情况，hP自适应法将在下文进行介绍。

在现有技术中已经实现了若干方法以更准确地并且更高效地对波传播系统进行建模。第一方法是以下方法，并且在图5和图6中示出。hP自适应FEM(以及hP自适应DGM)是采用具有基于后验误差估计器局部地调整的可变维度h和多项式次数P的元素的有限元方法(以及间断伽辽金方法)的普通版本。hP-FEM的由来追溯到Ivo Babuska等人的开创性工作[I.Babuska,B.Q.Guo:The h,p and h-p version of the finiteelement method:basis theory and applications,Advances in EngineeringSoftware,Volume15,Issue3-4,1992]。Babuska发现当使用h细化(将元素分为较小的元素)和P细化(增加元素的多项式次数)的合适的组合来细化计算方案时，有限元方法快速地指数收敛。该逻辑与标准FEM计算非常不同，因为不是以给定频率组装和计算单个解，组装并且求解了具有增加的h或P复杂度的解的序列。在每次计算之后，后验误差估计器分析该结果并且定义新的阶和网格细化直到判断解符合要求。[J.R.Stewart,T.J.R.Hughes:An a posteriori error estimator and hp-adaptive strategyfor finite element discretizations of the Helmholtz equation in exteriordomains,Finite Elements in Analysis and Design,Volume25,Issues1–2,1997]。

在大部分情况下，hP-FEM方法基于分层函数空间。在这些空间中，通过仅添加新的形函数来根据基础B(P)获得较高阶的基础B(P+1)。这对于P自适应有限元代码是必不可少的，因为当增加给定元素的内插的阶时，不必完全改变其形函数。图4公开了阶＝3的四边形上的分层基础的示例。第一行包括未修改的阶P＝1的形函数。

hP自适应方法在计算上非常高效，因为其使得能够真正地优化对于给定问题的计算工作量。然而，不足之处是需要自动网格细化。另外，在编程方面，hP-FEM是非常具有挑战性的。与标准的FEM解器相比，hP-FEM更难以实现和维持。

因此，一些现有技术方法提出简单地使用高阶有限元方法，而不采取任何自适应细化方案(不讨论先验估计器或后验估计器)，并且脱离多频率计算的背景。这通常被称为P-FEM方法，如在以下参考中的那样：[S.Petersen,D.Dreyer and O.von Estorff:Assessment of finite and spectralelement shape functions for efficient iterative simulations of interioracoustics,Computers Methods in Applied Mechanics and Engineering,195,6463-6478,2006]和[P.E.J.Vos,S.J.Sherwin and R.M.Kirby:From hto p Efficiently:Implementing finite and spectral/hp elementdiscretisations to achieve optimal performance at low and high orderapproximations.Journal of Computational Physics,229,2010]。

另一方法是以下方法。如已经解释的，在大多数的声学应用中，驻波传播问题的解需要在一个以上的频率上进行。往往需要跨越很宽的频率范围(音频频率范围)并且应对大量的频率。这类问题被称为频率扫描问题。加速该多频率问题的方法被称为频率快速扫描(FFS)法。FFS法试图加速标准低阶FEM而不是求助于较高阶的解。FFS法在几个具体的频率对解进行求解，并且基于如在以下参考中的解及其高阶导数的知识在周围近似地推断解[M.S.Lenzi,S.Lefteriu,H.Beriot,W.Desmet:A fast frequencysweep approach using Padéapproximations for solving Helmholtz finiteelement models,Journal of Sound and Vibration,Volume332,Issue8,2013]。

本发明的目的是通过在精度、计算时间和占用内存方面选择更高效的方法来优化在宽的频率范围内的系统的声音发射和传播的模拟。

发明内容

因此，在第一方面，根据权利要求1，本申请涉及用于预测在宽频率范围内的声源的发射以及其在周围介质内的传播的处理。

本发明立足于hP自适应方法和快速频率扫描方法二者的交汇。在本发明中，开发了新型的并且独创的方法以更高效地计算频率扫描问题，但是没有任何近似。本发明提出了一种算法，用于使用基于分层高阶离散化以及先验估计器的使用的域方法来解决在多个频率的波传播问题。

优点是，得益于高阶形函数的使用，工作量可适于每个频率的需要(在低频率的低工作量和在高频率的高工作量)。

又一个优点是，得益于先验估计器的使用，每个频率只计算一个解，而不是像现有技术的hP自适应FEM那样计算解的序列。

最后，得益于形函数基础的分层，在频率循环之前，能够预先计算与所有元素相关的质量基本矩阵K_e,i和刚度基本矩阵M_e,i，而不考虑内插的阶随频率改变的事实。

所有这些方面使得能够显著地降低多频率模拟所需的计算时间。

本发明的第二方面通过根据权利要求7的软件程序产品来实现。

当结合附图时，通过下文的详细描述，本发明的其他优点和新特征将变得更加清楚。

附图说明

根据下面的描述和附图将更好地理解本发明。

图1示出了汽车空腔内部的声音发射的示例，灰度等级表示在给定频率的空腔内部的可听得见的声压的水平；

图2示出了低频f₁所需要的汽车的稀疏网格；

图3示出了高频所需要的图2的汽车的精细网格；

图4示出了在上文和下文中提及的给定阶P的分层基础；

图5示出了关于给定的2D模拟问题的典型hP自适应FEM的中间解；

图6示出了图5的问题的典型hP自适应FEM的最终解；

图7示出了关于在频率f₁的简单2D网格由先验估计器给定的元素阶的重新分配，其中呈现了元素e₁和e₂；

图8示出了关于在频率的简单2D网格由先验估计器给定的元素阶的相似的重新分配，其中呈现了元素e₁和e₂；

图9、图10和图11示出了本发明的优选处理序列的流程图，图9示出了最大元素阶P_e,max的计算处理，图10示出了基于P_e,max的所有基本矩阵的计算处理，图11示出了多频率计算的最终阶段；以及

图12示出了分层嵌入基本矩阵的示意性表示。

具体实施方式

在此所描述的本发明提出了分层高阶域方法(FEM或DGM)的高效实现方案，用于通过基于专用先验误差估计器的高效多频率策略来更高效地解决频率扫描问题。在本发明的上下文中没有考虑P自适应或h自适应方案，网格是固定的，并且每个频率仅执行一次计算。

使用有限元空间的分层性质，通常仅在感兴趣的最大频率计算一次弱公式化。然后子矩阵简单地被提取，并且基于先验误差估计器在较低频率处被组装。

该方法在于结合三种概念：

-使用先验误差估计器，

-使用分层以构造质量和刚度矩阵，

-多频率计算。

1)所使用的第一概念是先验误差估计器。在hP-FEM方法中，先验误差估计器有时用于在hP细化处理中帮助设计最初网格。先验误差估计器表示在方程求解之前将向有限元网格内部的每个元素分配给定的内插阶的预处理技术。它是基于以下项输出在给定频率f_i的关于给定的1D、2D或3D元素e(直线、三角形、四边形、四面体、六面体、棱柱)所需的阶P_e,i的函数：

-特征元素长度h_e(例如最大边缘长度)、

-频率值f_i、

-特征局部介质性质(声速如果存在的话，平均流量大小)。

在第一近似中，下式可以应用于接近P_e,i：

$P_{e,i}=6\frac{h_e{f}_i}{c_{0e}}$

该式保证每个元素的自由度的密度等于6。

借助于示例，图7和图8示出了如何在实践中采用该式以定义示例结构(矩形网格)中的元素阶。图7和图8分别表示在频率f₁和f₂的网格中的元素阶P_e,i的重新分配的示例，f₂高于f₁。通过使用先验估计器已确定元素阶P_e,i。这两个图表示元素阶直接与元素大小和频率值成比例。如果声速不是频率相关的(大部分时间是这种情况)，还表示通常将在感兴趣的最高频率处找到最大阶。

为清楚起见，在图7和图8中两个元素被编号为e₁和e₂。

对于这两个元素，

P_1，1＝3(频率f₁下的e₁)

P_2，1＝2(频率f₁下的e₂)

P_1，2＝10(频率f₂下的e₁)

P_2，2＝9(频率f₂下的e₂)。

2)在本发明中所使用的第二概念是分层。通过低阶方法，相比于系统求解(全局系统矩阵的分解)，系统组装(基本矩阵的评估)相当廉价。通过高阶方法，这两种主要操作之间的效果平衡被修改，并且由于在基本方程中存在高阶积分，所以基本矩阵的评估可以变成与系统求解一样是计算密集的。

在此，没有使用分层，用于：

-后验增加关于给定元素的阶来改进在一个频率处的解(与hp自适应方法中的操作相同)，但是

-鉴于加速对多个频率的计算，从一个频率到另一频率先验增加整个网格的阶。

3)本发明中所使用的第三概念是多频率计算。其指代如下事实：在宽频率范围内进行该计算。

在图9、图10和图11中示出了用于在宽频率范围内预测声源的发射以及其在周围介质中的传播的处理的实施方式。

由用户提供的输入列表如下：

-适合于表示问题的有限元或间断伽辽金方法的N_e个元素的网格，

-定义问题的边界条件、源和/或材料性质的集合，

-升序的N_f个频率的列表。

该处理开始于可以被分成两个阶段的预处理部分：

-最大元素阶评定，

-最大基本矩阵计算。

预处理部分的第一阶段是最大元素阶评定。在该阶段期间，通过调用网格中的每个元素的先验估计器将元素阶P_e,i与每个元素e和每个频率f_i相关联。然后，在频率范围内针对每个元素确定最大元素阶P_e,max。这可以以数学形式被写为：

$=P_{e,\max }=\underset{1}{\max }{P}_{e,t}.$

返回我们关于图7和图8的先前示例，该最大元素阶评定的阶段将导致：

P_1，max＝10，以及

P_2，max＝9。

在图9中表示该第一阶段的完整处理。

在该处理的第一步骤(1)中，读取给定为输入的网格，并且对于每个元素获得特征元素维度h_e。

对于每个元素e和每个频率f_i：

-在步骤(2)中，引入局部流体性质。

-在步骤(3)和步骤(4)中，调用先验估计器以确定与每个元素相关联的最大元素阶P_e,max。在流体性质不是频率相关的情况下，将根据分析的最大频率处的值确定最大阶。在频率相关流体材料的情况下，可能出现其他场景。

在预处理部分的第一阶段结束时，对于每个元素以大小N_e的阵列来存储最大元素阶P_e,max。

预处理部分的第二阶段是最大基本矩阵计算。

在图10中呈现了该第二阶段的完整处理。

对于每个元素e：

-在步骤(5)中，加载最大阶元素P_e,max。

-然后，在步骤(6)中，通过使用最大元素阶P_e,max来分别形成频率无关的基本质量和刚度矩阵M_e,max和K_e,max。

-连续地，在步骤(7)中，以二进制格式将这些基本矩阵M_e,max和K_e,max存储在磁盘上。

重要的是注意，得益于先验估计器，针对每个元素，对于与最大阶元素相对应的阶，仅计算一次基本矩阵M_e,max和K_e,max。在该阶段没有频率循环。

在这两个预处理阶段完成之后，多频率计算部分可以开始于在用户所需的频率上的循环。在图11中呈现了该第二部分的完整处理。

对于每个频率f_i：

-对于每个元素e：

o首先，在步骤(8)中，计算该频率和该元素所需的元素阶P_e,i。

o然后，在步骤(9)中，调取与P_e,max相对应的最大基本矩阵K_e,max和M_e,max。

o在步骤(10)中，从K_e,max和M_e,max提取与内插阶P_e,i相对应的基本矩阵K_e,i和M_e,i。使用分层对该步骤至关重要。在现有技术中，在每个频率从无开始计算K_e,i和M_e,i矩阵，然而在本发明中，由于基础的分层结构，包含与P_e,i相对应的矩阵作为P_e,max的矩阵的子集。这在图12中示出，其中应当注意的是，较低阶的矩阵被包含在较高阶的矩阵中。

o在步骤(11)中，将基本矩阵K_e,i和M_e,i组装成全局矩阵K_i和M_i，它们分别表示系统在频率f_i的刚度和质量。在步骤(11)中，还存储活跃的自由度列表(该特定频率所需要的未知数的子集)。

然后，在关于元素的循环的实现之后，对于每个频率f_i：

-在步骤(12)中，与该特定频率f_i相对应的矩阵被组装

Z_i(f_i)＝K_i-(2πf_i)²M_i+C_i(f_i)

如果需要频率相关特征(如边界条件)，则将它们的贡献添加在C_i(f_i)中。通常Z_i在低频处很小(其中P_e,i很小)，并且当频率接近最大频率时较大。因此，对于每个频率该计算结果为隐含适应的。

-在步骤(13)中，表示源的右手侧被组装。

-在步骤(14)中，然后使用线性求解器来对该系统进行求解。

上面所描述的方法的优点基于如下事实：基于先验误差估计器在每个频率适应效果，并且基于如下事实：在频率循环之前仅对基本矩阵评估一次，这可观地减少了多频率模拟所需的计算时间。

Claims (7)

1.一种用于在宽频率范围内预测声源的发射以及声音在周围介质内的传播的处理，其中包括所述源和所述周围介质的系统由小且简单的元素e表示，

对于每个元素e和每个频率f_i，

-通过先验误差估计器使参数P_e,i与所述元素相关联(3)，所述参数P_e,i描述用于近似在频率f_i的声场的多项式次数的特征，

-确定关于所述元素的参数P_e,max(4)，其对应于所述估计器计算的在所述频率范围内的最大P_e,i参数，

对于每个元素e，

-使用所述参数P_e,max确定基本矩阵K_e,max和M_e,max(6)，它们描述所述元素分别对所述系统的刚度和质量的贡献的特征，

对于每个频率f_i，

o对于每个元素e，

■参数P_e,i用来确定用于近似所述声场的多项式次数(8)，

■对于所有元素，从所述矩阵K_e,max和M_e,max提取基本矩阵K_e,i和M_e,i，并且将所述基本矩阵K_e,i和M_e,i组装成分别表示所述系统的刚度和质量的全局矩阵K_i和M_i(10)(11)，

■基于全局矩阵K_i和M_i来建立全局矩阵系统Z_i(12)，

■使用线性求解器对所述全局矩阵系统进行求解(14)。

2.根据权利要求1所述的处理，其中，将表示所述系统的网格给定为所述处理开始时的输入。

3.根据权利要求1和2中任一项所述的处理，其中，将必须被采样的所述频率范围内的离散频率列表给定为所述处理开始时的输入。

4.根据权利要求1至3中任一项所述的处理，其中，将系统的边界条件、源和材料性质的集合给定为所述处理开始时的输入。

5.根据权利要求1至4中任一项所述的处理，其中，对于每个元素e引入局部流体性质(1)。

6.根据权利要求1至5中任一项所述的处理，其中，所述全局矩阵系统Z_i具有下面的形式K_i：

Z_i(f_i)＝K_i-(2πf_i)²M_i+C_i(f_i)

K_i和M_i分别表示所述系统的刚度和质量，C_i(f_i)是因所述边界条件出现的所有其他频率相关项，并且f_i是所关注的频率。

7.一种用于执行根据权利要求1至6中任一项所述的处理的软件步骤的计算机产品。