跨声速机翼阻力分解的一种改进方法

白策 包芸 张怀宝 王光学

摘要：提出黏性区域探测器的一种改进形式，并用于捕捉激波和翼梢涡的熵增阻力;给出尾迹平面的可压缩涡动力学诱导阻力表达式，并与基于热力学的诱导阻力对比。在跨声速来流状态下，对ONERA M6和某民用飞机巡航状态下的机翼阻力进行分解，同时分析该民用飞机机翼安装翼梢小翼前、后的远场阻力构成。结果表明：新的区域探测器合理可靠，黏性阻力与伪熵阻力的计算结果更加准确;2种诱导阻力计算方式的计算结果一致，但基于涡动力学的诱导阻力计算方法受积分平面位置的影响更小;安装翼梢小翼基本不影响整个流场的黏性阻力，减阻的主要效果体现为诱导阻力的减小。

關键词：跨声速机翼; 翼梢小翼; 尾迹积分; 诱导阻力; 熵增阻力

中图分类号：V211.41; TP211.3

文献标志码：B

An improved drag decomposition method for transonic wing

BAI Cea BAO Yuna ZHANG Huaibaob WANG Guangxueb

（a. School of Aeronautics and Astronautics; b. School of Physics Sun Yatsen University Guangzhou 510006 China）

Abstract：

An improved viscous zone detector is proposed to capture the entropy drag of shock wave and wingtip vortex. The induced drag calculation based on compressible vortex force on the wake plane is derived. The induced drag is compared with the thermodynamical one. The wing drags of ONERA M6 and a civil aircraft in cruise state are decomposed under transonic incoming flow. The farfield drag composition of the civil aircraft is analyzed before and after the winglet installation. The results show that the new zone detector is reasonable and reliable and the calculation results of viscous drag and spurious entropy drag are more accurate. Both of the two induced drag calculation methods can obtain consistent results and the induced drag based on vortex force is less affected by the integral plane position. The installation of the winglet does not affect the viscous drag of the entire flow field and the main effect of drag reduction is the reduction of induced drag.

Key words：

transonic wing;winglet; wake integral; induced drag; entropy drag

0 引 言

跨声速机翼设计特别重要，阻力因数的细微差别可能会对燃油消耗产生很大影响，增大运营成本和温室气体排放。MEREDITH[1]研究表明，对于典型的商用飞机，减少1个阻力因数单位（0.000 1）相当于多运载一名乘客（91 kg），这意味着更低的能耗和更优的经济效益，可以大大提高飞机的环保特性和市场竞争力。因此，对飞机在巡航阶段的飞行阻力构成进行精细化分析尤为重要。

在CFD计算中，采用压力和黏性剪切力在物面积分的方式，可以获得飞机的压差阻力和摩擦阻力，却无法从阻力的产生机制上进一步分析;远场阻力分解技术通过求解远场面积分或流场体积分的方式，可以从阻力来源上分析阻力构成。

传统远场法是由VAN DER VOOREN等[2]从动量守恒方程推导而来的，通过分析热力学可逆与不可逆过程，把总阻力分解成熵增阻力（主要是激波和黏性引起的）和诱导阻力。PAPARONE等[3]为Oswatitsch熵增公式增加2阶项，并分割流场域，将熵增阻力进一步分解为激波阻力、黏性阻力和伪熵阻力（数值误差）。在计算诱导阻力方面，BETZ[4]首次在测量二维阻力的风洞试验中提出尾迹积分法。不久，该理论被推广到三维风洞试验阻力测量中，其主要思想是在垂直来流的方向设置积分平面，通过积分尾迹面上无关熵增的物理量得到基于热力学的诱导阻力。WU等[5]和LIU等[6]从普朗特升力线理论出发，建立远场气动力与兰姆矢量积分的联系，从而得到基于涡动力学的诱导阻力。该方法现已推广至三维可压缩流动计算中。 [711]

阻力分解方法广泛应用于研究和工程中。LANZETTA等[12]研究黏性流动中超大远场对阻力分解精度的影响;YAMAZAKI等[13]将阻力分解技术应用在翼梢小翼多目标优化中;陈真利等[14]、高飞飞等[15]和李典等[16]采用阻力分解技术分析客机标模的阻力构成;VOS等[17]分析DPW4网格因素对阻力分解结果的影响。基于热力学和涡动力学的阻力分解技术已经扩展到非定常流动中。[1819]

本文选用“中场法”分解熵增阻力，给出基于尾迹平面的可压缩涡动力学诱导阻力表达式，并与基于热力学的诱导阻力进行对比。对ONERA M6机翼进行阻力分解，提出黏性区域探测器的一种改进形式。对某民用飞机跨音速巡航机翼安装翼梢小翼前、后的远场阻力构成进行分析，并比较2种诱导阻力计算方式的优缺点。

1 阻力分解理论基础

1.1 远场和近场阻力平衡理论

将x轴设置为来流方向，定常情况下对包含无动力物面的流场域进行面积分，得

式中：

式（2）左边为压力和黏性应力在物面上的面積分，对应近场法阻力Dnear的压差阻力Dp和摩擦阻力Df;等式右边为远场法远场阻力Dfar的计算式。将远场变量引入远场阻力以减小计算波动，定义

式中：下标∞表示流场变量在远场边界的取值。因此，远场阻力可化简为

根据热力学关系式可以导出

uu∞=fΔpp∞，ΔsR，ΔH

[WTHX]u[WTBX]2∞=

1+2ΔHu2∞-2γ-1Ma∞Δpp∞+1γ-1γeΔsR

γ-1γ-1

（5）

式中：γ为比热容;R为气体常数;Ma∞为来流马赫数。由此可见，速度受压力变化Δp、熵的变化Δs和焓的变化ΔH的影响。将流向速度变化分解成不可逆熵增过程和可逆过程2部分，其中不可逆熵增过程产生的阻力为熵增阻力，包括激波阻力Dw、黏性阻力Dv和伪熵阻力Dsp;可逆过程产生的阻力即诱导阻力Di。Dnear和Dfar的各阻力成分之和在理论上应该是相等的，即阻力平衡公式为

Dp+Df=Dw+Dv+Dsp+Di

（6）

1.2 熵增阻力的分解

对u/u∞应用泰勒级数进行线性化展开，可以得到熵增过程导致的速度变化为

Δ=f1ΔsR+f2ΔsR2+

OΔpp∞3，ΔsR3，ΔHu2∞3

（7）

式中： f1=-1γMa2∞; f2=-1+（γ-1）Ma2∞2γ2Ma4∞。

因此，不可逆熵增部分对阻力的贡献，即熵增阻力为

PAPARONE等[3]采用熵的差分平衡方程推导认为，熵增阻力的被积函数实际是熵增阻力的当地增长率。因此，理论上熵增阻力在物理熵变化为0的区域应该为0。然而，在CFD计算中不仅有物理熵增，还会因为数值算法的耗散和临近远场边界的网格粗化引入非物理熵增，必须分离这些区域。将整个积分区域

Ω拆分为Ω=Ωv+Ωw+Ωsp，Ωv为包含边界层和黏性尾迹的黏性区域，Ωw为包含激波的区域，Ωsp为剩余区域（定义为伪熵区域）。在各区域分别积分，即可得到熵增阻力的各成分，即

Dw=∫Ωw

[WTHX]Δ[WTBX]

·[WTHX]f[WTBX]irrevdΩ

Dv=∫Ωv

[WTHX]Δ[WTBX]

·[WTHX]f[WTBX]irrevdΩ

Dsp=∫Ωsp

[WTHX]Δ[WTBX]

·[WTHX]f[WTBX]irrevdΩ

（9）

分离各区域所用的探测器将在第2节的跨声速验证算例中详细介绍。

1.3 诱导阻力的计算

VAN DER VOOREN等[2]通过远场力的矢量形式给出诱导阻力的体积分形式。

由于黏性矢量与物面垂直，即

这样，可得到诱导阻力的体积分形式为

积分区域可取距物面距离为d的空间。该方法可较好地满足远、近场阻力平衡条件，但计算结果受积分区域影响，且难以应用于试验测量。在足够远的尾迹下游设置与来流方向垂直的积分平面，通过平面积分计算诱导阻力，更符合工程实际。[4]在积分平面内有

在线性远场中，边界法向量

上述方法基于远、近场阻力平衡，将速度变化进行热力学分解，得到热力学上的诱导阻力。还有一种方法是建立涡动力学的兰姆矢量与诱导阻力的联系。WU等[5]研究证明，在不可压缩流动中，旋涡力为扰动速度兰姆矢量的体积分，诱导阻力源自该旋涡力在来流方向的投影。随后，LIU等[6]为可压缩流动添加修

正项-

距离相关的矢量），从而将体积分限制在黏性和激波区域。MELE等[9]还指出，在马赫数大于0.7时，修正项

[WTHX]m[WTBX]ρ不再产生效果。为克服数值解在激波区域精度低的问题，将激波区域旋涡力的体积分等效替换成面积分

LIU等[20]进一步证明，在可压缩流动中，该旋涡力DLIU的形式为

式中：u′为扰动速度。

该面积分表达式证明，在足够远的线性远场内，诱导阻力项将退化消失。然而，尾迹区域的速度密度变化仍然明显，可在尾迹下游设置垂直于来流方向的平面，从而将基于涡动力学的诱导阻力DiV表示为积分平面面积分的形式，即

式中：

[WTHX]n[WTBX]∞为自由来流方向的单位向量。速度扰动为0的区域积分结果也是0，因此可将积分区域限制在与尾迹中心距离为d的区域内，但需要注意计算结果应对d的取值不敏感。

2 算例验证与区域探测器的选择

2.1 ONERA M6机翼的近场阻力求解

选取ONERA M6机翼外形进行跨声速算例验证，计算软件为结构网格求解软件SunFlow，选用SpalartAllmaras湍流模型。采用HO型结构网格划分计算域，物面网格见图1a）。远场距物面的距离约为120倍平均气动弦长，网格数量约为600万个，来流马赫数Ma∞=0.84，攻角α=3.06°，雷诺数Re=11.71×106，机翼65%展向站位的压力因数分布与试验数据的对比见图1b），其中x/clocal为无量纲弦长。计算结果与文献[21]试验数据基本吻合。

a）机翼表面和空间网格分布

b）65%展向站位压力因数分布

2.2 激波区域探测器

对于激波区域，采用基于压力梯度的探测器隔离激波的无黏部分。由黏性产生的熵比由激波产生的熵大几个数量级，因此即使黏性熵的微小污染也会大大扭曲激波阻力的计算结果。定义激波探测函数为

式中：Ma为当地马赫数。Fshock在可压缩区域满足Fshock>0，在激波区域满足Fshock>1，故Fshock可以捕捉激波上游区域。激波探测函数等价于马赫数在压力梯度方向的投影，使用RankineHugoniot關系式可以预估激波下游的马赫数。计算截止常数Kcw用于选择激波的上、下游区域，区域探测函数F满足Fshock>Kcw时认为是激波区域。Kcw的定义为

式中：Fshock max为Fshock的最大值。

LOVELY等[22]提出改进判据F′shock，即

F′shock增加额外项考虑激波的移动，适用范围更广泛。激波区域的探测结果见图2，其中ds/R为当地熵增。探测器捕捉到机翼的λ激波，且没有混杂边界层的熵增。

2.3 黏性区域探测器

定义边界层和黏性尾迹区域的区域探测函数为

式中：

μl和μt分别为层流黏性系数和湍流黏性系数。[3]当Fviscous满足Fviscous>Kbl·Fviscous，∞时，认为所在位置属于边界层和黏性尾迹区域，其中Kbl为截断常数（通常设为1.1），Fviscous，∞为探测函数在远场边界的取值。要注意区域的划分应对Kbl取值不敏感。

未加修正项的黏性探测器可视化结果见图3a）。在近似二维流动的翼根位置，采用基于Fviscous的探测器可较好地识别出边界层和尾迹区域，基本覆盖局部熵增较大的区域。然而，在靠近翼梢的区域，基于Fviscous的探测器只捕捉到边界层和尾迹部分的熵增，没能捕捉到翼梢涡对上翼面影响导致的熵增。MHEUT等[23]提出uirr表达式成立的条件为p

为正确捕捉物理熵增区域，对传统边界层和尾迹探测器添加修正项，该修正项选取当地总压损失为标准，仅在流动变化剧烈的区域起作用，其表达式为

Fviscous，fix=bfix·1-P0ρ∞a2∞

（21）

式中：bfix为用于调整修正项比重的常数。修正后的探测器将激波及波后的熵增区域一并纳入，因此实际捕捉的是物理熵增作用明显的区域。该区域减去

F′viscous得到的激波区域即为黏性区域，记新的黏性区域探测器为F′viscous，其可视化结果见图3b）。

a）未加修正项的探测结果

b）添加修正项并排除激波区域的探测结果

2.4 阻力分解结果随积分平面的变化

为减少网格数量，边界层和尾迹区域网格划分较密集，远场边界附近区域网格较稀疏。远场的数值耗散要大于实际情况，因此实际计算时通常只对积分平面的尾迹区域进行积分，以排除数值耗散和网格粗化的影响。

积分平面法得到的阻力构成随dx/c的变化曲线见图4。

图4中：DffT为基于热力学的远场总阻力;DffV为基于涡动力学的远场总阻力;Dv为积分平面上游的黏性阻力计算结果，用以直观表示尾迹区域黏性阻力的增加;下标old和new分别指Fviscous和F′viscous的熵增阻力分解结果;dx为积分平面距机翼后缘的距离;c为平均气动弦长。

由此可以看出，随着dx/c的增大，黏性阻力增大，诱导阻力减小，这

是诱导阻力在黏性尾迹区域的耗散的结果。使用Fviscous的确会将一部分黏性阻力当做伪熵阻力，导致黏性阻力项偏小，而改进的F′viscous可正确识别黏性区域，从而获得更准确的黏性阻力和更小的伪熵阻力。

远场阻力约在dx/c=60的位置与近场阻力平衡。诱导阻力的2种计算方法结果基本一致，只是DiV在dx/c<1时，出现诱导阻力计算偏小的情况，这是因为积分平面的流动变量不能满足线性化条件。另外，DiV在1

3 翼梢小翼效果分析

对某民用飞机机翼安装翼梢小翼，计算和分析巡航状态下的阻力构成。由于巡航阶段耗时最长，基于该状态的减阻设计对燃油经济性的提升最明显。该机翼的外形见图5a），网格划分方式与ONERA M6机翼算例保持一致，网格数量约为2 000万个。选用SpalartAllmaras湍流模型，在同升力条件下对原机翼和带翼梢小翼构型机翼进行数值模拟，取CL=0.5，来流马赫数Ma∞=0.785，攻角

αYJY=2.891 5°、αDXY=2.836 7°，雷诺数Re=25.3×106，远场距物面的距离约为90倍的平均气动弦长。SunFlow的近场计算结果表明，加装翼梢小翼后，原机翼阻力因数减小约29.9个阻力因数单位，近翼面的流线分布见图5b）。

a）机翼外形对比

b）上翼面边界层的流线对比

由此可以发现，翼梢小翼有整理翼梢气流的作用，有翼稍小翼机翼的翼梢涡气流的掺混程度明显小于原机翼。

2个构型机翼激波与黏性区域的探测结果见图6。在翼梢的尾迹位置，无翼梢小翼机翼翼梢涡的当地熵增更强，影响范围更大;安装翼梢小翼后翼梢涡当地熵增明显减小，尾迹区域的判定范围也小于无翼梢小翼机翼。

a）无翼梢小翼机翼的激波区域

b）无翼梢小翼机翼的黏性区域

c）有翼梢小翼机翼的激波区域

d）有翼梢小翼机翼的黏性区域

2个构型机翼的远场阻力构成随积分平面位置dx/c的变化见图7。与ONERA M6机翼算例一样，基于F′viscous的区域划分结果可大大提升黏性阻力和伪熵阻力的准确性。将体积分区域限制到dx/c=0的位置，记此时F′viscous得到的黏性阻力为Dv0，代表翼面边界层区域对阻力的贡献。原始机翼的Dv0约为67.0个阻力因数单位，安装翼梢小翼后，对应的Dv0有所增加，为87.6个阻力因数单位，但单独积分小翼部分的黏性阻力仅为3.2个阻力因数单位，说明Dv0的增加不仅是因为翼梢小翼

增加额外的面积，还与主翼面的流动状态发生改变有关。

a）无翼梢小翼机翼

b）有翼梢小翼机翼

将体积分区域拉伸至远场时，可用黏性阻力与Dv0的差值Ddv表示黏性尾迹区域的物理熵增对熵增阻力的贡献。原始机翼的Ddv为30.3个阻力因数单位，安装翼梢小翼后Ddv仅为14.8个阻力因数单位，可量化翼梢小翼对尾迹气流的整理效果。

综上所述，翼梢小翼有增加边界层熵增阻力、减小尾迹区域熵增阻力的效果。安装翼梢小翼机翼的黏性阻力仅增加约5.0个阻力因数单位。

计算诱导阻力时，随着积分平面向流场下游移动，部分诱导阻力会随着黏性尾迹的影响转化为熵增阻力。可在飞行器后缘附近设置积分平面计算机翼近场的诱导阻力，采用线性化插值方法近似得到貼近机翼后缘的诱导阻力：原机翼和有翼梢小翼机翼的热力学诱导阻力DiT分别为119.7和92.5个阻力因数单位，2个构型的涡动力学诱导阻力DiV分别为116.5和88.7个阻力因数单位。基于涡动力学的诱导阻力结果更有参考价值，这归功于其近场计算结果受积分平面的影响较小。

然而，将飞行器后缘附近求得的诱导阻力作为远场阻力分解的结果并不合适，因为飞行器的黏性尾迹一直延伸至流场下游远场，在计算黏性阻力时，在黏性尾迹的耗散部分已经包含诱导阻力，这导致黏性阻力和诱导阻力的计算很难匹配。由图7可以发现，对于不同机翼构型，很难选择相同的积分平面位置使得每个构型都能满足远场和近场阻力平衡条件，难以对不同构型的阻力成分进行对比分析。本文选取满足远场和近场阻力平衡条件较好的积分位置作为阻力分解的最终结果，黏性阻力和伪熵阻力选取整个流场基于改进探测器的区域划分结果，见表1。由此可知，安装翼梢小翼后原机翼的黏性阻力增加5.0个阻力因数单位，激波阻力减小1.9个阻力因数单位，减阻的主要效果体现为诱导阻力的减小，约减小30.0个阻力因数单位，占原机翼诱导阻力的29%。

4 结 论

将改进的远场阻力分解技术应用于结构网格求解软件SunFlow的计算结果后处理中，改进区域探测器，对黏性区域、激波区域和伪熵区域进行更精确

地划分。在诱导阻力计算中，将

传统的尾迹积分法引入

到基于可压缩涡动力学的诱导阻力计算中，并与基于热力学的诱导阻力进行对比。对某民用飞机机翼安装翼梢小翼前、后的阻力构成进行量化分析，研究各阻力成分随积分平面位置的变化情况。结论如下：

（1）改进的探测器对各区域的划分符合当地的实际熵增情况。本文提出的黏性探测器改进形式可以准确捕捉到激波后部和翼梢涡处的局部熵增，更适用于跨声速流场，可提高黏性阻力计算的准确度，减小伪熵阻力的计算结果。

（2）2种诱导阻力计算方法在远场积分平面均得到合理可靠的结果。基于涡动力学的诱导阻力计算结果受积分平面位置的影响较小，尤其是近场的积分结果一致性更好，因此考察飞行器近场诱导阻力时，基于涡动力学的诱导阻力更具有参考价值。

（3）翼梢小翼有整理小翼尾迹气流的效果，可减小黏性尾迹区域的熵增阻力，但增大边界层的熵增阻力，略微减小激波阻力。综合来看，翼梢小翼对整个流场熵增阻力的影响较小。翼梢小翼减阻的主要效果体现为诱导阻力的减小。

参考文献：

[1]

MEREDITH P T. Viscous phenomena affecting highlift systems and suggestions for future CFD development[C]// AGARD Conference Proceedings 515： HighLift System Aerodynamics. Alberta： Advisory Group for Aerospace Research & Development，1993.

[2] VAN DER VOOREN J SLOOFF J W. CFDbased drag prediction： State of the art theory prospects： NLR TP90247U[R].

[3] PAPARONE L TOGNACCINI R. Computational fluid dynamicsbased drag prediction and decomposition[J]. AIAA Journal 2003 41（9）： 16471657. DOI： 10.2514/2.7300.

[4] BETZ A. A method for direct determination of wingsection drag： NACATM337[R].

[5] WU J Z MA H Y ZHOU M D. Vorticity and vortex dynamics[M]. Heidelberg： Springer 2006： 2135. DOI： 10.1007/9783540290285.

[6] LIU L Q WU J Z SHI Y P et al. A dynamic counterpart of lamb vector in viscous compressible aerodynamics[J]. Fluid Dynamics Research 2014 46（6）： 061417. DOI： 10.1088/01695983/46/6/061417.

[7] MARONGIU C TOGNACCINI R. Farfield analysis of aerodynamic force by lamb vector integrals[J]. AIAA Journal 2010 48（11）： 25432555. DOI： 10.2514/1.J050326.

[8] MELE B TOGNACCINI R. Aerodynamic force by Lamb vector integrals in compressible flow[J]. Physics of Fluids 2014 26（5）： 545556. DOI： 10.1063/1.4875015.

[9] MELE B OSTIERI M TOGNACCINI R. A novel midfield breakdown of aerodynamic force in compressible flows[C]// Proceedings of 34th Applied Aerodynamics Conference. Naples： AIAA 2015. DOI： 10.2514/6.20163429.

[10] MELE B OSTIERI M TOGNACCINI R. Vorticity based breakdown of aerodynamic force in threedimensional compressible flows[J]. AIAA Journal 2016 54（4）： 11981208. DOI： 10.2514/1.J054363.

[11] YAMAZAKI W MATSUSHIMA K NAKAHASHI K. Aerodynamic design optimization using dragdecomposition method[J]. AIAA Journal 2008 46（5）： 10961106. DOI： 10.2514/1.30342.

[12] LANZETTA M MELE B TOGNACCINI R. Advances in aerodynamic drag extraction by farfield methods[J]. Journal of Aircraft 2015 52（6）： 18731886. DOI： 10.2514/6.20123211.

[13] YAMAZAKI W MATSUSHIMA K NAKAHASHI K. Drag prediction decomposition and visualization in unstructured mesh CFD solver of TAScode[J]. International Journal for Numerical Methods in Fluids 2008 57（4）： 417436. DOI： 10.1002/fld.1643.

[14] 陳真利 张彬乾. 基于尾迹积分的阻力计算方法研究[J]. 空气动力学学报 2009 27（3）： 329334. DOI： 10.3969/j.issn.02581825.2009.03.012.

[15] 高飞飞 郝海兵 颜洪. 阻力精确分解的方法研究与数值模拟[J]. 航空计算技术 2017 47（2）： 2932. DOI： 10.3969/j.issn.1671654X.2017.02.008.

[16] 李典 郝海兵 颜洪 等. 大型客机远场阻力分解方法验证研究[J]. 航空计算技术 2018 48（5）： 510. DOI： 10.3969/j.issn.1671654X.2018.05.003

[17] VOS J SANCHI S GEHRI A. DPW4 results using different grids including nearfield/farfield drag analysis[C]// Proceedings of 28th AIAA Applied Aerodynamics Conference. Illinois： AIAA 2010. DOI： 10.2514/6.20104552.

[18] GARIEPY M TREPANIER J Y MALOUIN B. Generalization of farfield drag decomposition method to unsteady flows[J]. AIAA Journal 2013 51（6）： 13091319. DOI： 10.2514/1.J051609.

[19] OSTIERI M TOGNACCINI R BAILLY D et al. Aerodynamic force and lamb vector field in compressible unsteady flows[C]// Proceedings of 2018 AIAA Aerospace Sciences Meeting. Florida： AIAA 2018. DOI： 10.2514/6.20180548.

[20] LIU L Q SU W D KANG L L et al. Lift and drag in threedimensional steady viscous and compressible flow[J]. Physics of Fluids 2017 29（11）： 126. DOI： 10.1063/1.4989747.

[21] SCHMITT V CHARPIN F. Pressure distributions on ONERAM6 wing at transonic mach numbers： AGARD AR138[R].

[22] LOVELY D HAIMES R. Shock detection from computational fluid dynamics results[C]// Proceedings of 14th Computational Fluid Dynamics Conference. Norfolk： AIAA 1999. DOI： 10.2514/6.19993285.

[23] M HEUT M BAILLY D. Dragbreakdown methods from wake measurements[J]. AIAA Journal 2008 46（4）： 847862. DOI： 10.2514/1.29051.