一种非紧致可渗透边界气动噪声数值积分计算方法与流程

本发明属于流体力学和声学技术领域，涉及一种非紧致可渗透边界气动噪声数值积分计算方法。

背景技术：

近年来，航空飞行器、民航飞机、高速列车、叶轮机械噪声问题大量出现，因而对气动噪声数值研究的需求不断增大。复杂结构表面往往存在尖角或者内凹外凸等构型，如三段翼模型、风机叶片模型，不易开展气动噪声数值研究工作。实际应用中，往往采用流动高精度计算，所带来的计算量增加了计算负担。然而，低精度流动计算又无法获得准确的声源信息。尤其边界形状较为复杂时，传统计算手段存在数值奇异性问题，无疑增大了数值计算的难度。

如边界的几何尺寸远小于声波波长，可将边界看作是紧致的，声波传播到边界时将继续向前辐射，可看作自由空间内的声传播；反之，如边界的几何尺寸大于或接近于声波波长，则边界可看作是非紧致的，声波向空间传播时，会产生向外辐射的声波，同时在边界上发生散射，产生二次辐射声波。非紧致结构的声传播与自由空间空间的声传播存在着明显的差异。

边界积分方程方法结合边界元的思想，有效提高了气动噪声问题的计算效率。近年来，大量学者采用该类方法对非紧致结构开展了诸多研究工作，基本分为两种，一种是寻求满足给定边界条件的精确格林函数解，该方法通过构造格林函数的方法获得格林函数并结合低精度流动计算进行声学计算；第二种是结合压强变量分解表达式，提出远场声压计算的声学积分方程，通过捕获物面边界上的散射声源进而完成远场噪声计算。

当前相关研究工作，对于简单边界较为适用，复杂边界情形下存在计算效率低、计算复杂度高的问题。

技术实现要素：

本发明的目的是在于克服现有技术计算方法中，采用物体表面为积分边界，结合边界元方法实施积分计算的过程中存在计算效率低、计算复杂度高、数值奇异性的缺点与不足，开展复杂结构气动噪声数值计算的研究工作，提供一种非紧致可渗透边界气动噪声数值积分计算方法，该方法采用包围物体的可渗透边界为积分边界，可以有效提高数值计算效率、降低计算复杂度并克服数值奇异性问题。

本发明上述目的通过以下技术方案予以实现：

一种非紧致可渗透边界气动噪声数值积分计算方法，

包括以下步骤：

步骤1：采用流动计算软件(如fluent)将流场计算区域离散为k个网格单元，获得声源信息，包含密度ρ、压强ph、速度ui(i＝1，2，3)，对于低马赫数流动，ρh≈ρ；

步骤2：在物体附近区域选择边界sp，并将sp离散为l个网格单元，一般选择为与流动离散网格相同的边界；

步骤3：采用公式计算sp上等效的散射声源pa(z，ω)，可重新整理为线性方程组：其中l＝1，2，....l，e是单位矩阵，h为对角线为零的对称矩阵

式中，m＝1，2，3，…l，n＝1，2，3，…l，l为边界s包含的单元数目；

步骤4：采用公式计算远场监测点x的声压pa(x，ω)。

上述计算表达式中，ω代表圆周频率，c0代表声波传播速度，vk代表第k个流场网格单元的面积或体积，[]k代表第k个流场网格单元上的信息，[]l代表第l个边界单元上的信息，pa(zl，ω)代表边界上点zl在ω频率的声压，pa(x，ω)代表远场监测点x在ω频率的声压，最后一项中j代表虚数单位。δij代表二阶张量，

进一步，当所述问题为二维模型，步骤3、4实施过程中，上述二维模型通过下列公式计算自由空间格林函数g(ym，yn，ω)

进一步，当所述问题为三维模型，步骤3、4实施过程中，通过下列公式计算自由空间格林函数g(ym，yn，ω)

确定ym点和yn点的自由空间格林函数沿yn点的外法向偏导数包括：

当所述问题为二维模型，步骤3、4实施过程中，通过下列公式确定

当所述问题为三维模型，步骤3、4实施过程中，通过下列公式确定

确定ym点和yn点的自由空间格林函数沿yn点的二阶偏导数

包括：

当所述问题为二维模型，步骤3、4实施过程中，通过下列公式确定

进一步，当所述问题为三维模型，步骤3、4实施过程中，通过下列公式确定

确定x点和yn点的自由空间格林函数沿yn点的外法向偏导数包括：

当所述问题为二维模型，步骤3、4实施过程中，通过下列公式确定

当所述问题为三维模型，步骤3、4实施过程中，通过下列公式确定

确定x点和yn点的自由空间格林函数沿yn点的二阶偏导数

包括：

当所述问题为二维模型，步骤3、4实施过程中，通过下列公式确定

进一步，当所述问题为三维模型，步骤3、4实施过程中，通过下列公式确定

其中，下标m，n＝1，2，3，…，m为网格单元编号，为第一类0阶hankel函数，第一类1阶hankel函数，k＝ω/c0为声波数。

通过步骤3计算物面散射声源，在右端项求解完成的前提下，可以通过计算

获得矩阵

h的子项hmn，待求声源点ym所对应的矩阵列项为

其中s＝s1∪s1…∪sl，ym点与yn重合时，hmn＝0。

通过步骤4计算远场监测点声压，在右端项求解完成的前提下，可直接获得。

本发明的优点及有益效果：

1.本发明采用二阶精度流动计算获得流场声源信息，有效提高了数值计算的工作量；

2.本发明通过计算可渗透边界上的声源获得等效散射声源，降低了复杂物面边界散射声源的计算复杂度问题；

2.本发明通过计算辐射噪声和散射噪声获得远场噪声，声压信息中可区分辐射部分和散射部分；

3.本发明可以解决复杂构型气动噪声或多个物体流场内的气动噪声传播计算，适用性较强。

附图说明

图1是流场计算局部网格图；

图2是点源声波传播示意图；

图3是可渗透边界及物体边界位置分布图；

图4是层流圆柱流场计算区域；

图5是层流圆柱圆柱附近网格密度分布；

图6是低频率下层流圆柱边界积分计算的远场声压指向性，其中，(a)f＝f0，(b)f＝2f0；

图7是层流圆柱表面及可渗透边界的位置分布；

图8是不同积分边界获得的远场声压指向性，其中，(a)f＝f0，(b)f＝2f0；

图9是湍流圆柱表面及可渗透边界的位置分布；

图10是不同积分边界获得的远场声压级分布，其中，(a)f＝f0，(b)f＝2f0。

具体实施方式

下面结合附图和具体实施方式对本发明的技术方案作进一步详细地说明。

本发明中涉及的声场计算建立在流动计算的基础上，以二维空间内圆柱和翼型的声学计算为例，如附图1所示，将每一个网格看做一个声源，远场检测器接收到的声音信号为所有声源发出的声音信号之和。

已知点声源传播遵循波动方程(1)，

每个流场声源可以等价看作一个点源，同时遵循流动传播规律，满足lighthill波动方程(2)

考虑到流动的可压缩性，将密度ρ，压力p和速度u分解为如下形式

通过一系列积分推导，可以获得如下频域空间内声学计算模型

其中，ω＝2πf，表示圆周频率；z点为监测器位置，光滑积分边界对应的c(z)＝1/2，远场监测点对应的c(z)＝1，v0表示物面边界与积分边界之间的流场区域，v-v0表示积分边界外部的流场区域。

由上述声学计算模型可知：

1.在流动计算完成基础上，需沿时间方向存储流动计算获得的每个网格单元密度ρh，压力ph和不同坐标方向的速度ui，uj。

2.为了监测声学信息，监测器可放置在任意位置，一般情形下包含两类：物体表面和非物体表面，其中非物体表面可以为流动计算区域内或者流动计算区域外，

以单个声源的传播为例，如附图2所示，在a点放置声源点，在c点放置监测器，b为物体表面边界，c点接收到的声音信号包含两部分：(1)a→c；(2)a→b→c

其中，由途径(1)获得的噪声为辐射声，由途径(2)获得的噪声为散射声。流动计算完成的前提下，可直接计算辐射噪声；物面b分布的散射声源计算完成后，才可计算散射噪声。对于实际流动问题，由模型方程出发，选择s为物体表面边界可开展声学计算。

然而，当物体表面结构较复杂或者流场内存在多个物体时，以二维问题为例，如图3中包含翼型和圆柱两个对象，翼型边界和圆柱边界分别为sc和sd。采用物体表面边界进行积分计算时，积分边界s＝sc+sd，散射声源计算需要将物面边界离散为细密的边界单元，存在计算复杂度高、数值奇异性等问题。考虑到上述问题，采用包围物体的可渗透边界sp进行积分计算，此时模型方程(4)中s＝sp，具体实施步骤为：

1：采用流动计算软件(如fluent)将流场计算区域离散为k个网格单元，存储所有网格单元随时间变化的流场信息，即声源信息；

2：在物体附近区域选择边界sp，并将sp离散为l个边界单元，一般选择为与流动离散网格相同的边界，以减少计算工作量；

3：采用公式(5)计算sp上等效的散射声源pa(z，ω)，可重新整理为线性方程组(6)。

4：采用公式(7)计算远场监测点x的声压pa(x，ω)。

公式(5)如下所示

公式(6)如下所示

其中l＝1，2，....l，e是单位矩阵，h为对角线为零的对称矩阵

式中，m＝1，2，3，…l，n＝1，2，3，…l，l为边界s包含的单元数目。

公式(7)如下所示

当所述问题为二维模型，通过下列公式计算自由空间格林函数g(ym，yn，ω)

当所述问题为三维模型，通过下列公式计算自由空间格林函数g(ym，yn，ω)

确定ym点和yn点的自由空间格林函数沿yn点的外法向偏导数包括：

当所述问题为二维模型，通过下列公式确定

当所述问题为三维模型，通过下列公式确定

确定ym点和yn点的自由空间格林函数沿yn点的二阶偏导数

包括：

当所述问题为二维模型，通过下列公式确定

当所述问题为三维模型，通过下列公式确定

确定x点和yn点的自由空间格林函数沿yn点的外法向偏导数包括：

当所述问题为二维模型，通过下列公式确定

当所述问题为三维模型，通过下列公式确定

确定x点和yn点的自由空间格林函数沿yn点的二阶偏导数

包括：

当所述问题为二维模型，通过下列公式确定

当所述问题为三维模型，通过下列公式确定

其中，下标m，n＝1，2，3，…，m为网格单元编号，为第一类0阶hankel函数，第一类1阶hankel函数，k＝ω/c0为声波数，c0为声波在介质中传播的速度。

采用公式(6)计算物面散射声源，在右端项求解完成的前提下，可以通过计算

获得矩阵

h的子项hmn，待求声源点ym所对应的矩阵列项为

其中s＝s1∪s1…∪sl，ym点与yn重合时，hmn＝0。

采用公式(7)计算远场监测点声压，在右端项求解完成的前提下，可直接获得pa(x，ω)。

进一步，分析声源各项的贡献：

采用公式(6)进行边界散射声源计算，右端第一项

表示计算单元的声源信号的辐射贡献，右端第二项

表示流场区域v-v0内所流场网格单元声源信号对z的辐射贡献，右端第三项

表示边界s上其它单元的声源信号对z的辐射贡献，最后一项

表示流场区域v0内声源信号对z的等效辐射声，从而获得所有声源对边界观测点的贡献，即散射声源。

采用公式(7)进行远场监测点声压时，右端第一项

表示流场区域v-v0内所有流场网格单元声源信号对x的辐射贡献，右端第二项

表示边界s上声源信号对x的辐射贡献，右端第三项

表示边界s上散射声源对x的贡献，即为散射噪声。最后一项

表示流场区域v0内的声源信号对x的等效辐射声，从而获得所有声源对远场观测点的贡献。

通过具体应用证明本发明有以下特点：

本发明采用包围物体的可渗透边界为积分边界，通过计算可渗透边界上的声源获得等效散射声源；远场噪声中可以有效区分辐射噪声和散射噪声；

本发明突破传统方法采用物面边界计算散射声源的界限，可以获得近场等效散射声源；同时可以考察近场流动区域的噪声贡献；

本发明可以解决复杂构型气动噪声或多个物体的流场内气动噪声的传播计算，适用性较强。

本发明无需采用高精度计算，大大提高了计算效率。

本发明选择包围所有物体的光滑流场网格边界为可渗透边界，容易实施。

具体应用案例：

实施例1层流圆柱噪声预测

选取直径d＝1m的二维圆柱为研究对象，来流马赫数ma＝0.15、雷诺数re＝100。本小节主要考察非紧致可渗透积分边界的可靠性。流场计算区域及网格分布如附图4所示，圆柱中心与坐标原点重合，流场计算区域限定在水平方向[-12d，24d]和竖直方向[-12d，12d]范围内，共包含44，080个四边形网格单元，u代表来流速度，u＝ma×c0。局部区域网格分布图见附图5。为与khalighi的dns方法计算结果进行数值比较，以(1.86d，0)为圆心、r＝12.9d为半径作圆周，在圆周上均匀布置180个观察点。

首先，为验证计算方法的有效性。采用圆柱边界(cylinderwall)计算散射声源。附图6给出了低频率下的声压指向性图，其中红色圆圈代表khalighi采用dns(直接数值模拟)的计算结果，蓝色实线代表当前计算方法采用圆柱表面进行积分计算获得的结果，两种方法吻合。

为了验证可渗透边界的有效性，选择如附图7所示的可渗透边界s1、s2、s3为积分边界，附图8给出了不同积分边界计算等效散射声源获得的远场声压指向性。由附图8可知，不同可渗透边界进行积分计算获得的计算结果与圆柱表面相当吻合。可渗透边界积分计算方法可以精确的捕捉等效的散射声源。其中，图8中，f0代表漩涡脱落频率。

实施例2湍流圆柱噪声预测

本算例主要分析湍流噪声分布，以三维圆柱模型为研究对象，物理参数选取如下：圆柱直径d＝0.019m，来流马赫数ma＝0.2，对应re数为9×10⁴。以圆柱中心为坐标原点，计算分析x向(-15d，20d)、y向(-15d，15d)和z向(0，4d)范围内的流体运动诱发的噪声。附图9给出了圆柱表面及不同可渗透边界的位置分布图。

附图10给出了不同积分边界获得的远场声压级分布。由图10可知，当f＝f0、f＝2f0，声波的最大辐射方向分别表现为竖直与水平方向，这与hu等基于格林函数数值解的边界元方法获得的结果一致，也符合在漩涡脱落频率(f＝f0)和二次谐波频率(f＝2f0)气动噪声主要呈现偶极子分布的物理事实。三个积分边界都能够准确的漩涡脱落频率(f＝f0)的噪声分布。不同之处在于，当选择s3为积分边界，在二次谐波频率(f＝2f0)，左侧声压级偏小，说明实际工程问题中需要选择合适的可渗透边界。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。