一种适用计算声学问题的网格间断人工黏性分布策略

王 鑫， 余培汛， 杨 海， 潘 凯

(中国飞机强度研究所 航空声学与动强度航空科技重点实验室， 西安 710065)

抑制短波的有效手段主要有两种：① 迎风格式，通过固有的耗散性来抑制短波，该方法的不足之处在于无法识别短波长波范围，将两者一起耗散，这样影响到计算结果的准确度；② 在对称模板的差分方程中添加人工黏性项，所选择的人工黏性项必须在有效抑制短波的同时，保证对长波的影响很小。Tam在1993年提出了人工选择耗散法[3]，通过对人工黏性耗散项的系数选择，使得其耗散项只对短波起作用，而对长波几乎没有影响，用于气动声学计算时取得了很好的结果。但是该方法只限于均匀对称网格，尚未在非均匀网格中进行拓展。Engquist[4]提出用非线性滤波法来捕获激波的方法，它在每一时间步用一个非线性过滤因子对数值解进行处理,只作用在部分网格点上。王泽晖等[5]在Bjrn提出的非线性滤波法基础上，根据一个检验函数，分析了非线性滤波效率与波长的关系，并引入了具有TVD性质的非线性滤波法，通过数值算例验证了DRP格式与非线性滤波法结合的可行性，表明该方法对于在网格交界面、初始值间断、非线性声波传播等计算气动声学常见情况下产生的短波有很好的抑制作用。该方法对短波具有很好的抑制作用，但同时对长波的耗散也相对较大，影响数值的计算精度。

本文在直角坐标人工选择耗散法的基础上，将其转换成曲线坐标系下的表达形式，并在多块网格交界面处、交界点、出口边界等处按照当地网格雷诺数的倒数进行正态分布，并通过NASA标模算例进行了该方法的验证分析。

1 数值模拟方法

1.1 线化欧拉方程

本文的主动控制方程采用的是二维线化欧拉方程作为声传播方程，其表达形式如下

(1)

式中:Q0=(ρ′,u′,υ′,p′)T为瞬时变化量(密度，速度分量和压力)。Q0=(ρ0,u0,υ0,p0)T为时均变量。式(1)的矩阵系数表示为

(2)

式中:γ=1.4为比热比，S为声源项。

经过曲线坐标转换后，式(1)可写为下式

(3)

1.2 空间离散格式

计算气动声学所要模拟的是物理波的传播过程，为此所选择的差分格式要能准确反映其传播特征(色散性，耗散性，方向性，相速度，群速度等)。高阶泰勒级数截断法所获得的差分模板关注的是声波的耗散性。它能准确的逼近物理波的幅值，保证波幅计算值与实际值偏差较小。但对于声传播问题，差分格式不仅需要关注声波的幅值，而且需要关注如何准确描述声波的传播，降低振荡，避免发散。

本文在对于式(1)中通量项的离散采用了Tam提出的7点DRP差分格式，其x方向的通量项可表示为式(4)，其余两个方向的通量采用同样的差分格式。

(4)

其中,插值系数为

(5)

图1 有效波数随波数变化曲线(Tam方法)

1.3 无反射边界条件

在计算气动声学中，由于受到计算资源的约束，计算区域不可能无限大，因此需要人为地在计算区域中采用一种边界条件来取代远场区域的信息。这种取代远场区域的边界条件，其需要起到两个作用：① 让计算区域内的各种物理波能无反射的通过边界；② 让边界以外的任何扰动不能传入计算区域，防止污染数值解。

论文中采用了Hu所提出的PML(Perfectly Matched Layers)边界条件，对无分裂形式的PML方程[7-9]在曲线坐标上进行了转化，并改善边界条件的数值离散格式，增强其稳定性。

远场边界处的二维LEE(Linearized Euler Equations)方程组通过拟线性化处理，可写为

(6)

根据Hu的方法，可知由于对流效应，波的相速度和群速度可能沿着不同的方向发展。如果此时在PML方程中仅引入数值黏性，反射波将不能完全被消除，以致计算的稳定性受到影响。为了解决这一难题，Hu采用了保持空间不变而对时间进行变换的方法，即：

(7)

(8)

将其代入到式(6)中，可得以下方程

(9)

(10)

然后通过复数变易法，引入阻尼，即：

(11)

将式(11)代入式(10)可得

(12)

式(12)两端同时乘以(ω+iσξ)(ω+iση)/ω2，可得式(13)，即：

(13)

傅里叶逆变换得

(14)

(15)

将上述偏导数关系式(15)代入式(14)，展开可得：

(σξ+ση)Q+σξσηq+βQ(σξξxA1+ηxA2ση)+

β(ξxA1+ηxA2)σξσηq=0

(16)

2 曲线坐标系下人工黏性在间断处的分部策略

对于声学问题的数值模拟，在划分网格的过程中不可避免地带来网格间断等问题。而这些间断(物面边界、远场边界、网格交界面)易引起差分格式的色散，产生数值伪波。以物面边界为例，对于空间离散采用7点DRP格式，假如物面不采用虚网格的形式，如图2所示。垂直于物面的网格点A、B空间离散将无法采用对称模板的7点DRP格式。为此在添加人工黏性项时，网格点A、B的人工黏性模板也是有所区别的，对于物面上方网格点A采用3点人工黏性模板，网格点B采用5点人工黏性模板，网格点C和D采用7点模板。网格交界面处可通过采用3层虚网格添加7点模板的人工黏性项，也可类似于物面边界方式处理(注：程序在网格交界面采用了3层虚网格的形式)。

另外为了保证计算的稳定性，在公共边的奇异点(非常规的1-1交接网格顶点处)附近需要添加额外的人工黏性项,如图3所示。接下来将介绍程序中在网格内部、网格交界面、物面边界、公共边奇异点处所采用的人工黏性项分布方式。

图2 间断处的人工黏性函数分布示意图

图3 奇异点示意图

假设采用7点人工黏性模板，在一维控制方程中，添加人工黏性项后，可表示为

(17)

其中RΔx=a0Δx/υa，为单位网格雷诺数。当网格存在拉伸和变形时，即计算坐标与直角坐标不一致，存在坐标变换时，式(18)中人工黏性的添加需要转换到计算坐标系下

(18)

如果保证式(17)与式(18)中人工黏性添加的效果相同，需使得RΔξ=xξRΔx。Δξxξ为曲线坐标中长度Δx在计算坐标中相对应的值。与一维方法类似，在三维空间中添加人工黏性可表示为

(19)

(20)

式中:ξ0为网格边界处坐标;nΔξ在文中取为当前网格块在ξ方向的半宽。同理在其余两个方向网格雷诺数的倒数分布表示为

(21)

(22)

在需要添加人工黏性的奇异点(ξ0,η0,ζ0)上网格雷诺数的倒数分布为

(23)

3 算例验证

为了验证曲线坐标系下人工黏性格式的合理性，文中选用了NASA第二届气动噪声会议[11]中初始扰动绕圆柱传播算例。同时，该算例同时也验证了远场无反射边界条件处理曲面边界的能力。

如图4所示，直径r=1的圆柱的坐标系原点位于圆柱的中心，点A、B、C分别是三个声压观测点，坐标分别为(r=5,θ=90°)、(r=5,θ=135°)、(r=5,θ=180°)。在距离原点r=4处分布了一个人工点声源S，其初始赋值如下

(24)

图4 圆柱计算模型示意图

圆柱壁面采用的是滑移边界条件，远场采用的是无反射边界条件。对于控制方程的离散，采用7点模板的DRP格式，人工黏性的分布策略如图5～图7所示，从网格块边界向内部网格雷诺数的倒数逐渐降低。文中计算网格的最小尺度为0.01，声场计算采用统一的无量纲时间步长Δt=0.002。图8(a)、图8(b)、图8(c)、图8(d)分别为t=50Δt、t=500Δt、t=2 000Δt、t=5 000Δt四个时刻的声传播布云图。从图8可知，声波随着时间推移，向空间逐步散开。当声波到达圆柱时，声波会发生部分反射现象，形成二次声源向空间传播。

图5 ξ方向无量纲分布

图6 η方向无量纲分布

图7 奇异点处无量纲分布

(a) t=50Δt

(b) t=500Δt

(c) t=2 000Δt

(d) t=5 000Δt

图8 不同时刻的声传播云图

Fig.8 Schematic diagram of sound propagation at different times

为了定量上对比分析所采用的人工黏性项分布方式的合理性，图9(a)、图9(b)、图9(c)分别表示观测点 A、B、C处的声压随时间的变化曲线。图9同时给出了计算解和解析解，其中圆圈表示解析解，实线为计算结果。由图9可知，文中计算得到的结果与解析解吻合得很好，对比结果说明了声波传播程序添加人工黏性项的合理性及处理曲边边界能力。

(a) A处

(b) B处

(c) C处

图9 观测点处声压随时间变化曲线

Fig.9 Curve of sound pressure at observation point

4 结 论

通过文中理论的阐述及算例的验证分析，可得出以下结论：

(1) 本文基于笛卡尔坐标系的人工黏性添加策略推导了曲线坐标系下的表示形式，并通过在网格间断等处的正态分布策略，提出了一种适用于网格间断处的人工黏性添加策略，通过标模算例的计算分析验证了该方法可成功应用于曲线坐标系下复杂网格的计算声学短波发散问题中。

(2) 针对初始扰动绕圆柱传播算例中的计算结果，可知声波首先向空间传播。当声波到达圆柱时，形成二次声源向空间传播。这一声传播现象以及监测点的对比分析，在验证网格添加策略合理性的同时也验证了远场PML边界条件的无反射特性。