旋翼干扰流场数值模拟中的新型Blend低速预处理方法

陈龙，夏健，田书玲

南京航空航天大学 航空宇航学院，南京 210016

旋翼/机身干扰和直升机/舰船干扰[1-2]等旋翼干扰非定常流场非常复杂。直升机机身和舰船流场会对旋翼气动力和操稳特性产生影响，旋翼尾迹对机身气动特性和舰面流场也有较大影响。流场中舰船及机身附近速度很低，是不可压流动。旋翼桨尖附近却处于跨声速状态并可能存在激波，是可压缩流动。因此在研究旋翼干扰流场时，需要考虑可压/不可压非定常流动的统一处理，应采用基于可压缩Navier-Stokes方程的低速预处理方法。旋翼干扰非定常流场的模拟若不采用低速预处理，过大的刚性和耗散会带来严重的收敛性和准确性的问题。低速预处理方法的核心问题是预处理矩阵的设计[3-6]。预处理矩阵对预处理方法的收敛性影响非常大，理想的预处理矩阵应具有特征值相对集中、消除刚性等特点。一般定常雷诺平均Navier-Stokes方程的预处理采用定常预处理矩阵[7-8]，非定常雷诺平均Navier-Stokes(URANS)方程的预处理采用非定常预处理矩阵[9-10]。

如果采用定常低速预处理模拟非定常问题，在双时间步长内迭代的收敛性问题会很严重[11]，甚至需要比无低速预处理更多的迭代。非定常低速预处理是在定常低速预处理的基础上引入矩阵中的伪时间和物理时间项，消除了刚性，在双时间步长内迭代中具有较好的收敛性。但是，非定常预处理在进行旋翼/机身干扰和直升机/舰船干扰为代表的旋翼干扰流场模拟时，物理时间步长需要按照旋翼转速确定。该时间步长较机身或舰船特征时间小很多，导致非定常低速预处理矩阵中的非定常截断参数大于1，这时非定常低速预处理自动关闭[12]。也就是说，在处理此类问题时，非定常低速预处理会自动退化到无预处理形式，收敛性和预测精度均无法得到保证。采用定常和非定常低速预处理在旋翼干扰流场数值模拟时遇到了两难困境。鉴于此，Potsdam等[12]将非定常低速预处理和定常预处理特征值矩阵进行组合，建立了Blend低速预处理方法。该方法解决了非定常低速预处理会自动退化到无预处理的问题，减小了流场中不可压区域的耗散。Sachdev和Hosangadi[13]针对特征值矩阵中定常预处理部分采用对应的定常预处理矩阵，改进Potsdam的方法为精确的代数形式，并推广到AUSM (Advection Upstream Splitting Method)类的通量分裂格式中。

本文在Potsdam和Sachdev的Blend低速预处理方法基础上结合对流项和黏性项量级分析，将定常低速预处理和非定常低速预处理特征值矩阵进行组合，建立出一种新型Blend低速预处理方法。本方法仍采用非定常预处理矩阵，对定常预处理和非定常预处理的通量格式耗散项进行组合得到新的通量格式耗散部分。在南京航空航天大学自主开发的OVERU软件[14]中应用该方法，经验证本文发展的方法相比传统方法计算效率大幅提高，也具有较好的预测精度，适用于旋翼干扰非定常流场的数值模拟，工程应用前景广泛。

1 控制方程与数值方法

控制方程为任意拉格朗日欧拉(ALE)形式的URANS：

(1)

(2)

式中：ρ为密度；u、v、w为速度；E为内能；W为守恒变量；xj为坐标；t为时间；Fj和Gj为通量。

采用非结构混合网格有限体积方法求解控制方程。空间离散采用格点格式，可处理多种单元类型。湍流模型可采用Negative-SA一方程模型和Menter的剪切应力输运(SST)两方程模型。时间离散采用隐式LU-SGS方法。非定常计算采用双时间步长方法进行二阶精度的时间离散。采用隐式重叠网格方法[14]处理网格间的相对运动。

2 低速预处理方法

预处理矩阵的建立首先需要将守恒变量W形式的方程转换成原始变量Q的形式，即

(3)

(4)

(5)

将原方程中的雅可比矩阵替换为预处理矩阵Γ，则预处理Navier-Stokes方程可写为

(6)

本文采用Weiss和Smith的预处理矩阵[15]：

(7)

(8)

式中：Map为定常截断参数；c为声速。

预处理参数采用当地马赫数配合全局截断的方法：

(9)

(10)

式中：Mai为当地马赫数；Ma∞为来流马赫数。

为了和预处理之前的方程保持相同形式，以便于程序编写，采用守恒形式的预处理方程为

(11)

称ΓM-1为守恒形式的预处理矩阵。定常预处理矩阵可以直接应用到双时间步长的计算中[10]，但由于其矩阵特征值的病态，造成了双时间步长内迭代收敛性很差。因此，需要建立非定常预处理矩阵，引入伪时间，对于双时间步长方法有

(12)

对该方程物理时间项采用二阶向后欧拉格式离散，伪时间项采用一阶离散得到离散方程，t是物理时间，τ是伪时间，有

(13)

式中：n表示第n个物理时间步；m表示第m个伪时间步。方程右端项为

(14)

令

(15)

则ΓuM-1为守恒形式的非定常预处理矩阵[11]。

令

(16)

(17)

Θ′=bΘ-(b-1)ρp

(18)

非定常预处理参数的定义为

(19)

(20)

式中：Mau为非定常截断参数；L为参考长度；Δt为时间步长。当高Strouhal数非定常流模拟时，时间步长Δt非常小；当模型的参考长度非常大时，将造成Mau>1，从而导致Map=1，此时非定常低速预处理为关闭状态。也就是说在处理此类问题时，非定常低速预处理会自动退化到无预处理形式，收敛性和预测精度均无法得到保证。旋翼机身干扰流场就属于此种情况。同时，对于复杂算例，特征长度和特征时间很难准确定义，这种不确定性影响着Mau的取值，也影响着非定常低速预处理的预测精度和收敛性。

3 Blend低速预处理方法

为解决非定常截断参数Mau>1造成非定常预处理自动关闭的问题，Potsdam组合定常低速预处理和非定常低速预处理特征值矩阵，从JST(Jameson-Schmidt-Turkel)格式的预处理形式出发建立了Blend低速预处理方法。

文献[12]认为不可压流场定常计算时，无预处理黏性项的速度场和温度场耗散过大，压强场耗散不足。通过定常低速预处理使得黏性项速度场、温度场和压强场均为合理的耗散。非定常计算时，问题就变得复杂的多，尤其是较小的物理时间步长或高Strouhal数的问题。若采用定常低速预处理模拟非定常问题，会造成病态特征值和双时间步长内很差的收敛性，同时人工黏性项中压强场耗散也会过大。若采用非定常低速预处理，时间步长较小造成Mau>1，非定常低速预处理就会自动退化到没有预处理的形式。为摆脱这个两难困境，Potsdam通过将定常低速预处理和非定常低速预处理特征值矩阵组合建立了Blend低速预处理方法[12]，方法中特征值矩阵为

(21)

(22)

引入选择矩阵

(23)

式(22)可表示为

(24)

在此基础上，Sachdev针对特征值矩阵中定常预处理部分采用对应的定常预处理矩阵改进了Postsdam的方法(式(25)方框项)，这种形式更适合通量差分格式耗散项的构造。文献[13]经过验证后证明Sachdev的方法和Potsdam的方法计算结果差别极小。Sachdev的Blend低速预处理中JST格式的人工黏性项为

D1/2=

(25)

量级分析方法是低速预处理收敛性研究的一个重要手段[3,16]。文献[3]从数值分析的角度，进行了对流项和人工黏性项的量级分析，研究了定常预处理方法对低马赫数时Euler方程收敛性的影响，当Ma→0时，未采用预处理的对流项的量级为

(26)

黏性项的量级为

(27)

当Ma→0时，预处理后对流项的量级为

(28)

黏性项的量级为

(29)

可见当Ma→0时，无预处理的对流项和人工黏性项量级并不匹配，有且仅有速度场的黏性项偏大，而在预处理后各方程量级保持一致。采用文献[3，16]的量级准则对Potsdam的Blend方法进行量级分析，在Mau>1时预处理系统的特征值量级为

(30)

(31)

则本文的Blend低速预处理中选择矩阵为

(32)

本文的Blend低速预处理中JST格式的人工黏性项为

D1/2=

(33)

在非定常预处理的框架下，当Mau>1时非定常预处理关闭，会出现格式耗散偏大的问题。而此时本文方法中对流项和黏性项的量级仍然是匹配的，减小了数值耗散，从而解决非定常低速预处理自动退化为无预处理的问题。同时本文的思想也适用于AUSM类格式的改进[13]。

4 计算结果及分析

4.1 圆柱非定常扰流

本算例模拟圆柱非定常扰流，取很小物理时间步长模拟非定常低速预处理自动退化为无预处理的问题，对本文建立方法的收敛性和正确性进行考核。网格采用三角形和四边形混合网格，约6万网格点，如图1所示。计算条件：来流马赫数为0.02，雷诺数ReD=1 000，采用层流计算，圆柱直径D=0.01 m，时间步长Δt=5×10-6s，非定常截断参数Mau=D/(πΔtc)≈1.87。此时非定常预处理自动退化为无预处理的形式。

图2给出了双时间步长内迭代的残差收敛曲线，可见定常预处理难以收敛，Blend预处理收敛性较好，本文方法略快于Sachdev的Blend方法，无/非定常预处理则因其耗散最大所以收敛速度最快。此处“无/非定常预处理”表示此状态Mau>1非定常预处理会自动退化到无预处理的形式，无预处理和采用非定常预处理计算结果是一致的。由于同样采用了非定常预处理矩阵，本文方法在Mau>1时预处理同样会关闭，但是改进过的黏性项与对流项量级仍然匹配，数值耗散也更合理。采样点(距圆心3.4D)的速度历程如图3所示，定常预处理、本文Blend预处理和Sachdev的方法几乎重合，而无/非定常预处理的过大耗散使其结果与其他计算结果差别很大。图2和图3中也对比了无/非定常预处理时采用Roe和HLLC (Harten-Lax-van Leer-Contact)迎风格式，可见无预处理时两种格式在此算例中无法收敛。图4和图5分别给出了无/非定常预处理和本文Blend预处理的涡量云图，可见无/非定常预处理的过大耗散使流场中的涡强在脱离圆柱后快速降低，而Blend预处理的涡量云图则显示出Blend低速预处理的耗散合理，得到了较好的涡量结果。综上，两种Blend预处理计算结果一致，但收敛效率较定常预处理有较大提高。而无/非定常预处理会耗散过大，计算结果较差。

图1 圆柱非定常扰流计算网格Fig.1 Computational grid for unsteady flow around cylinder

图2 圆柱扰流非定常残差收敛曲线Fig.2 Convergence history of residual of unsteady flow around cylinder

图3 采样点速度历程Fig.3 Velocity history of sampling point

图4 圆柱扰流无/非定常低速预处理计算结果Fig.4 No/unsteady low speed preconditioning result of flow around cylinder

图5 圆柱扰流本文Blend低速预处理计算结果Fig.5 Low speed preconditioning result of flow around cylinder of the proposed Blend method

4.2 旋翼/机身干扰流场

本算例进行直升机前飞时旋翼机身干扰的数值模拟，对本文Blend预处理方法收敛性和正确性进行考核。采用直径2.5 m的单旋翼直升机基准模型，包括4片桨叶的旋翼模型和Robin机身模型。旋翼实度为0.077，翼型采用NACA0012，弦长0.072 m，几何扭转0°，旋翼的周期变距为θ0=3.1°，θ1c=-0.6°，θ1s=1.5°，桨盘倾角αs=-5°，前进比为0.1。

图6为本算例所采用的重叠网格系统，共包含7块子网格，其中2块为远场背景网格，它们与绕Robin直升机机身的网格为静止网格，2块背景网格单元数分别为31万和163万，Robin直升机机身网格单元数为677万，均为六面体网格；绕每片旋翼叶片的网格随旋翼运动，网格单元数为107万。重叠网格系统总网格数约1 300万。物面第1层网格y+≤1。

图7给出了Robin直升机机身表面网格。湍流模型采用Negative-SA一方程模型。时间步长取旋翼转过0.5°方位角，参考长度为机身长度L=2.5 m，则非定常截断参数Mau=L/(πΔtc)≈38.8。此时非定常预处理自动退化为无预处理的形式。

图6 前飞直升机流场模拟重叠网格系统Fig.6 Overset grids system for flowfield simulation of helicopter in forward flight

图7 机身表面网格Fig.7 Fuselage surface grid

图8 定常预处理参数对单独机身气动力的影响Fig.8 Isolated fuselage forces vs steady preconditioning parameters

图9 单独机身压力系数分布计算值与试验值对比Fig.9 Comparison between experimental and computed distribution of isolated fuselage pressure coefficient

图10 单独机身压力云图对比Fig.10 Comparison of pressure contour of isolated fuselage

图11 定常预处理参数对收敛曲线的影响Fig.11 Convergence history vs steady preconditioning parameter

采用无/非定常低速预处理、定常预处理和本文Blend预处理方法开展直升机前飞时旋翼/机身干扰流场的模拟。随旋翼方位角的变化机身法向力系数的变化可以反映旋翼对机身气动力的干扰，图12和图13分别为采用不同双时间步长内迭代数和通量格式对机身法向力系数CN的影响。从收敛性方面考虑：无/非定常预处理时，机身法向力在内迭代(Nin为迭代次数)10、20和50次时很接近，仅需要约10步内迭代即可收敛；定常预处理时，从20增加到100次内迭代后机身法向力仍然没有收敛，收敛性非常差；本文Blend预处理，20次和50次内迭代的机身法向力几乎重合，表明其仅需约20次内迭代即可收敛，收敛性较定常预处理大幅改善。从流场预测精度考虑：充分收敛后定常低速预处理和Blend低速预处理计算结果表现出较好的一致性，无/非定常低速预处理由于存在较大耗散，JST格式和HLLC格式均得到较平缓的计算结果，存在较大的误差。

图14给出了不同方法的内迭代残差收敛曲线，可见本文Blend预处理方法的收敛效率略高于Sachdev的Blend预处理方法，定常预处理方法后期收敛速度会变慢，无/非定常预处理时JST格式由于其耗散最大收敛速度最快，HLLC格式则收敛困难。图15给出了旋翼方位角0°时各方法计算的机身表面压力分布，定常预处理方法和本文Blend预处理方法结果符合得很好，而无/非定常预处理结果中的压力分布存在振荡和较大误差。图16给出了流场的Q值等值面，本文Blend方法可以更好地进行旋翼尾迹的捕捉，而无/非定常预处理方法由于较大的耗散，显著影响了尾迹捕捉的精度。

图12 内迭代数对机身法向力的影响Fig.12 Fuselage normal force history vs number of sub-iterations

图13 通量格式对机身法向力系数的影响Fig.13 Fuselage normal force coefficient vs flux schemes

综上，本文建立的Blend低速预处理方法保持了定常低速预处理的预测精度，并大幅改善了收敛性，总体计算效率可提高约5倍，相比定常预处理方法可提高内迭代收敛标准，使计算精度得到保证。

图14 直升机前飞模拟时内迭代残差收敛曲线Fig.14 Residual value cnvergence history in simulation of helicopter in forward flight

图15 直升机前飞压力系数分布对比Fig.15 Comparison of pressure coefficients distribution of helicopter in forward flight

图16 直升机前飞流场的Q等值面Fig.16 Q criterion iso-surface for helicopter in forward flight

4.3 直升机/舰船干扰流场

本算例进行直升机/舰船干扰非定常流场的数值模拟，对本文Blend预处理方法收敛性和正确性进行考核。采用美国两栖攻击舰LHA的简化模型作为舰船模型，将4.2节直升机模型放大到18 m作为模拟的旋翼和机身模型，如图17所示。直升机位于舰首的第一起降点，桨盘中心距甲板13 m。舰船前进速度为20 m/s。旋翼的周期变距为θ0=10.3°，θ1c=-2.7°，θ1s=2.4°。

图18为本算例所采用的重叠网格系统，共包含6块子网格，分别为静止的舰船网格和直升机机身的网格，网格单元数分别为1 140万和677万。绕每片旋翼叶片的网格随旋翼运动，网格单元数为107万。舰船网格采用非结构混合网格，其余为六面体网格。重叠网格系统总网格数约2 440万 。物面第1层网格y+≤1。图19给出了直升机和舰船的表面网格。湍流模型采用Menter的SST两方程模型。时间步长取旋翼转过0.5° 方位角，参考长度为舰船甲板长度L=200 m，则非定常截断参数Mau=L/(πΔtc)≈489.9。此时非定常预处理自动退化为无预处理的形式。

图17 直升机/舰船干扰流场模拟模型Fig.17 Simulation model for helicopter/ship interaction flowfield

图18 直升机/舰船干扰流场模拟重叠网格系统Fig.18 Overset grids system for helicopter/ship interaction flowfield simulation

图19 直升机/舰船干扰流场表面网格Fig.19 Surface grids for helicopter/ship interaction flowfield

图20 直升机/舰船干扰模拟流程Fig.20 Flowchart for simulation of helicopter/ship interaction

图21 直升机/舰船干扰流场Q等值面Fig.21 Q criterion iso-surface for helicopter/ship interaction flowfield

图22 直升机/舰船干扰流场垂向速度Fig.22 Vertical velocity of helicopter/ship interaction flowfield

图23 直升机/舰船干扰模拟时残差收敛曲线Fig.23 Residual value convergence history for simulation of helicopter/ship interaction

5 结 论

本文针对定常和非定常低速预处理在旋翼干扰流场数值模拟时遇到的两难困境，在Potsdam和Sachdev的Blend低速预处理方法基础上结合对流项和人工黏性项量级分析，将非定常低速预处理和定常低速预处理特征值矩阵进行组合，建立出一种新型Blend低速预处理方法。

1) Blend低速预处理方法通过非定常低速预处理矩阵和定常低速预处理矩阵进行组合，解决了在旋翼干扰流场数值模拟时非定常低速预处理自动退化为无预处理的问题。在双时间步长内迭代中相比定常低速预处理收敛效率大幅提高，可提高内迭代收敛标准，使计算精度得到保证。

2) 无/非定常低速预处理在处理此类流场时，存在耗散过大的问题。充分收敛后定常低速预处理和Blend低速预处理计算结果表现出较好的一致性。

3) 本文建立的Blend低速预处理方法相比于Potsdam和Sachdev的方法，对流项和人工黏性项量级更匹配，修正了温度场黏性过小的问题，提高了收敛效率和鲁棒性。