在现如今经济快速发展的大环境下, 燃油价格仍然是高居不下, 而对于一架航班来说, 超过20%的运营成本都是消耗在了燃油使用上, 也因此整个航空工业界都在积极寻求创新的飞机设计方法来有效提升燃油使用效率^[1]。而且, 飞机油耗的减小和全机总阻力的降低息息相关, 即使全机阻力很小的提升也能为整个飞机有效载荷的增加提供重要的帮助。也因此诸如流动控制、气动优化和新型布局的探索等设计方法已经成为有效改善飞机气动性能的最主要手段, 而准确可靠的阻力预测方法是其必不可少的技术实现途径^[2]。

此外, 在当下的飞机设计中, 工程师也越来越依赖CFD方法进行阻力分析, 这是因为基于RANS方程求解的数值方法成功快速的发展和高可用度的计算资源实现, 使得现在针对复杂的全机构型进行流动计算和分析变得愈加高效方便。高精度的气动力数据对于飞机的精细化设计是必不可少的, 然而其仍旧依赖于湍流模型、转捩模型、网格质量和拓扑构造等^[3]。在这种程度下, 对于阻力的判定和量化就是非常重要的一件事。而根据工业界标准, 计算流体力学对于阻力系数预测的精度需要控制在1count以内。

另一方面, 通过应用空气动力学技术委员会(applied aerodynamics technical committee)的基金支持, NASA从2001年开始组织阻力预测大会(drag prediction workshop, DPW), 主要目的就是评估当下气动力预测的数值计算方法和工具。传统计算阻力的方法, 是通过积分机体表面的压力和剪应力来实现的, 通常称之为近场法。即使流场内局部求解是精确的(比如对于压力和速度型来说), 但对于计算结果来说仍然是不足以充分保证精确度的。尤其对基于Euler/RANS的数值求解而言, 通常会涉及到人工数值耗散和数值误差的出现, 从而会产生伪阻力(spurious drag)^[3], 这部分阻力的贡献在现有的计算资源和网格条件下很难可以做到忽略不计。另一个问题是, 近场法计算出来的阻力仅仅可以分解为压差阻力和摩擦阻力, 并不能将阻力按照物理形成原因进行分解, 比如黏性阻力、激波阻力和诱导阻力^[4]。

而另外一种方法, 是基于动量定理进行阻力预测的计算, 通常称之为远场法, 是一种可以解决上述问题的阻力预测方法。此方法是将阻力计算和熵的变化相结合。最早的先驱性研究工作是由van der Vooren和Slooff开展的, 然后由van der Vooren和Destarac逐渐发展成熟的。远场法对于气动工程师来说是一项实用功能很强的技术, 因为它不仅能够从阻力形成的物理机制上将阻力进行分解, 并可以将每部分阻力区域可视化, 而且能够量化伪阻力大小并将其从数值计算总阻力中移除, 从而有效提升阻力系数的计算精确度。

基于远场法计算阻力早在1925年由Betz^[5]首次提出, 而这个想法是来源于Von Kármán等^[6]给出的将守恒定律应用于有限控制体的思想。van der Vooren等^[7]提出了重要的阻力分解概念, 并针对远场法进行了系统综合的理论推导。van der Vooren等^[8]推导的基于热力学分解的阻力预测方法被广泛应用至今, 他们提出的方法不仅能够比近场法更精确地进行阻力预测, 还能将阻力按照形成的物理机制分解成激波阻力, 黏性阻力和诱导阻力, 此外还能计算得到因为数值耗散引起的伪阻力分量。而且, Destarac^[9]对于亚声速自由来流进行基于RANS方程的求解, 针对无动力飞机进行讨论(理论本身可以考虑带动力构型), 总结对比了不同公式下的计算结果精度, 并将阻力分解方法应用在不同的流场构型分析中。Paparone等^[4]提出了一种含有二阶修正项的方法, 推导出一个通用的远场法计算阻力的表达式。这种方法的基本思路和Destarac的方法是一致的, 只是他在求解轴向速度损失量时, 是基于泰勒级数展开的, 而且形成了一个关于熵, 压力和总焓等变量的函数。Gariépy等^{[1, 10-11]}的基础理论是在一个固定控制体中基于线性动量守恒关系式而形成的, 主要基于van der Vooren和Destarac的原始理论方法以及Laurendeau和Trépanier等人提出的改进方法, 并针对轴向速度损失量做出了相应的改进提升。此外, Marongiu等^[12]和Povitsky等^[13]开展了另外一派远场法研究理论, 他们用速度项代替压力项, 提出了涡矢量和Lamb矢量的概念。此方法主要是为了更好地表达诱导阻力分量, 但存在的问题是并不能将阻力清晰地分解为黏性阻力和激波阻力。

而国内相关学者也开展过一些方法研究:陈真利等^[14]利用尾迹积分法按阻力产生的物理机制, 将阻力分解为熵增阻力和诱导阻力进行数值计算。刘杰等^[15]采用尾迹面积分法, 通过在飞行器尾迹面上对相对流动参数进行积分得到升力和阻力, 并基于欧拉方程计算无黏流动下的升力、诱导阻力和激波阻力。李杰等^[16]采用基于动量定理的远场方法将飞机的阻力分解为黏性阻力, 激波阻力和诱导阻力, 并基于DLR-F4翼身组合体, 利用Euler方程求解对方法进行了验证。

基于国内外发展现状, 本文的研究目的和工作开展分为以下4个方面:①对比和分析近场法和远场法进行阻力预测的特点, 并提炼出现有主流的几种远场法关于轴向速度损失量的优劣势和差异, 提出本文关于轴向速度损失量的改进方法, 建立改进的基于远场法的阻力预测方法和阻力分解方法; ②针对阻力分解方法建立过程中, 需要对阻力区域的选择划分进行判断和决定的问题, 开展相关参数敏感性的讨论分析; ③基于本文改进的阻力分解方法, 针对CRM翼身组合体构型开展气动特性分析研究; ④将改进的阻力分解方法融入基于梯度的气动外形优化设计系统, 针对CRM构型开展气动外形优化设计研究。

1 阻力预测和分解方法 1.1 近场法阻力预测方法

在经典的近场法中, 阻力的表达式如下

(1)

式中：Body指机体表面；n是指向单元体外的单位正矢量；n_x=n·i且i是自由来流方向的单位矢量；p是压强, τ_x是黏性应力矢量。

(1) 式积分中的第一项和第二项分别对应压差阻力和摩擦阻力, 可见近场法的阻力分解很直接, 而且对于流动没有任何的假设。

1.2 远场法阻力预测方法

连续方程和x方向的动量方程表达形式如下

(2)

(3)

(2) 式和(3)式整合可得

(4)

此时定义矢量

(5)

设围绕机体的控制体V有边界∂V, 且∂V=S_A+S_F+S_D, 则有

(6)

在机体表面, 上式积分第一项为0, 因此有

(7)

(7) 式第一项即为近场法下阻力计算的表达, 因此(7)式可整理为

(8)

在边界S_F上, 流动可以认为是无黏无扰动的, 即τ_x=0, u=u_∞, p=p_∞, 且如果假设下游区边界S_D足够远, 则

(9)

此时, 定义轴向速度损失量

(10)

国内外现有主流的几种远场法计算阻力基本都是基于此理论体系下, 但是定义轴向速度损失量会导致各自理论方法有一定的优劣势和差异。

Destarac等^[9]对轴向速度损失量定义为

(11)

而Paparone等^[4]提出了一种含有二阶修正项的方法, 它基于泰勒级数展开, 是一个关于熵, 压力和总焓等变量的函数

(12)

式中, , , f_H1=1, , , , , , 。

此外, Gariépy等^[1]对轴向速度损失量定义为

(13)

1.3 改进的远场法阻力预测方法

综合1.2节几种不同的远场法, 可以得出其共同特点就是可以将热力学不可逆过程产生的阻力写成统一的形式, 即

(14)

而在Destarac的方法公式中对于Δu的表达, 即(11)式, 如果根式下为负值, 则积分就会出现一定问题, 从而导致阻力预测的不准确。

但是这种情况对于Paparone采取的泰勒级数展开就没有影响, 但是其方法又会存在截断误差的问题, 对于三阶项及以上采取了忽略的形式。

而Gariépy则是针对Destarac的方法存在的问题提出了自己的改进方法, 就是在根式下为负值时, 令Δu为(13)式的形式, 但是依然存在根式下需要判断正负的问题。

因此本文综合上述各种方法的优劣势和差异, 进行了新的补充和改善, 具体方法就是:首先采用Destarac的方法, 但是要针对网格单元进行(11)式中根号项内正负的判断, 正的区域就用其方法进行积分; 然后, (11)式中根号项内为负的区域采用Gariépy的方法, 但是要继续针对网格单元进行(13)式中根号项内正负的判断, 正的区域就用其方法进行积分; 最后, 剩下的区域采用Paparone的方法进行积分。

此时, 分别定义黏性阻力, 激波阻力和伪阻力

(15)

(16)

(17)

而诱导阻力的积分则延续Destarac的方法

(18)

(19)

此外, 本文追加定义远场法下的阻力分解量, 称为黏性压差阻力, 它是黏性阻力D_v=D_vp+D_f中由非摩擦阻力引起的阻力分量(由于流动分离和尾迹区等产生的), 因此也就有如下关系式

(20)

2 阻力分解方法的参数敏感性分析

在远场法预测阻力过程中, 针对不同的公式和实现方法, 有一些参数的选取会影响到最终阻力的预测和分解的精确度。本小节将会对这些参数的选取方法进行说明, 进而探索参数选取对最终阻力预测的敏感性影响。

2.1 区域划分算法与判断函数

从第1节对熵增阻力的划分可以看出, 阻力分解的准确性依赖于对黏性区域和激波区域的正确划分, 需要制定针对区域选择的算法规则。对于基于RANS的数值求解来说, 问题就是判断一个网格单元中的流体热力学信息是属于边界层尾迹区还是激波区, 亦或是空间中剩下的区域。

而在激波附面层干扰区域, 往往会存在区域判断而产生的分歧, 因此, 本文给定的判断算法规定:首先, 优先进行激波区域的判断; 然后, 再在余下的区域中进行黏性区域的判定, 最后, 剩下的区域都属于伪阻力积分区域。

2.1.1 激波区域判定

激波区域的判断识别方法依靠的是基于下面的无量纲判断函数

(21)

式中：a表示声速；▽p表示压力梯度。通常来说, 激波区域的选取需要依赖于阈值K_cw的大小, 整个流动区域满足F_shock>K_cw的地方即被选择为用来积分激波阻力的区域。

2.1.2 边界层尾迹区区域判定

本文选取了一种基于湍流流动的识别函数, 由于涡黏性系数是一种可靠的对黏性效应的反映, 所以选定的识别方法为

(22)

式中：μ_l和μ_t分别是层流黏性系数和涡黏性系数；F_bl的值在边界层尾迹区区域会非常高, 而在剩下区域F_bl≈1；黏性区域的选择通过F_bl>K_bl·F_bl∞来确定, 其中F_bl∞是自由来流的识别函数值, K_bl是对于黏性区域选取的阈值大小。

2.2 波阻区域的扩增和阈值判定

虽然激波判断函数和黏性区域判断函数能够起到识别作用, 但是他们有时并不能够很完美地选择所有应该识别的区域^[17]。对于激波阻力, 在很多情况下, 激波识别函数并不能准确地捕捉所有激波附近的波动, 这是由求解器的数值格式决定的。因此需要进行区域扩增来包含所有的激波阻力区域, 扩增的策略为:首先紧挨选择区域的单独一层网格被选作扩增区域的第一层, 然后剩下的几层会一层层紧挨着往外扩展, 等扩展到一定层数的时候, 如图 2所示, 之前未被捕捉到的激波阻力区域基本都会被涵盖。虽然这种方法看起来略显过度增加选择区域, 但是这些区域会继续加以判断函数的识别, 因此是合理可行的。

为了确定本文最终所采用的激波阻力区域扩展层数, 选用NACA0012翼型进行分析验证。验证求解的网格选取了3套计算网格, 分别是网格量为129×129的粗网格、257×257的中等网格和513×513的密网格。选择基于欧拉方程的数值求解方法, 计算状态为Ma=0.80, α=0°。图 3给出了激波阻力随扩增层数的变化, 从图中可看出, 随着扩展层数的增长, 激波阻力系数的预测结果趋向一致, 当选择扩展5层就可取得较好的结果。因此最终确定选取激波阻力区域的扩增层数为5。

此外, 激波识别函数是用来选取进行激波阻力积分的区域。当被判断网格区域的识别函数值大于阈值时, 这些网格就被选中用来积分^{[4, 17]}。而为了决定合适的识别函数阈值, 就需要分析阈值改变导致的阻力变化情况。

选用NACA0012翼型进行分析验证。选取粗网格进行验证求解。计算求解基于欧拉方程, 计算状态为Ma=0.80, α=0°。图 4给出了激波阻力数值随阈值K_cw的具体变化, 从图中可以看出, 随着阈值的增加, 激波阻力系数的预测结果变化较小, 当取K_cw就可取得较好的结果。因此最终决定在本文中选取K_cw=1.0。

2.3 黏性阻力区域的扩增和阈值判定

涡黏性系数在临近物面及附近区域值都是很小的, 所以这些地方会因为阈值的存在而容易未被捕捉到边界层区域中。因此需要针对物面附近的区域进行扩增选择^[17]。本文对黏性阻力区域进行扩增的过程如图 5所示。

本文选用RAE2822翼型进行分析验证。求解的网格选取了3套计算网格, 分别是网格量为25 600的粗网格、102 400的中等网格和409 600的密网格。选择基于RANS方程的求解器, 计算状态为Ma=0.734, C_L=0.824, Re=6.5×10⁶。图 6给出了黏性阻力随扩增层数的变化, 从图中可看出, 随着扩展层数的增加, 黏性阻力系数的预测结果趋于稳定, 当选择扩展4层就可取得较好的结果。因此最终选定黏性阻力区域的扩增层数为4。

此外, 黏性阻力识别函数也需针对阈值进行选取和分析^{[4, 17]}。仍选用RAE2822翼型进行了计算验证。求解的网格选取网格量为25 600的粗网格。计算基于RANS方程, 计算状态为Ma=0.734, C_L=0.824, Re=6.5×10⁶。从图 7中可看出黏性阻力和激波阻力随阈值的变化, 随着阈值的增加, 激波阻力系数和黏性阻力系数的预测结果均变化较小, 取K_bl=1.1就可取得较好的结果。因此决定在本文中选取K_bl=1.1。

2.4 诱导阻力的积分面选取

根据(18)式对诱导阻力定义可以看出, 诱导阻力积分依赖于对飞机下游区边界S_D的选取, 因此需要分析和确定诱导阻力的积分面选取。为此, 本文选用椭圆平面形状机翼构型进行计算验证, 其基本几何信息为:平面形状采用椭圆公式生成, 参考面积为3.06 m²。验证求解的网格选取了2套计算网格, 分别是网格量为192 512的粗网格和1 540 096的密网格, 图 14为构型的平面形状及粗网格分布示意图。而为了观察诱导阻力的变化情况, 本文选择欧拉方程在亚声速工况下进行求解, 计算状态为Ma=0.50, α=0°。

诱导阻力产生于流向动能转向横向动能的转换。远离物面时, 由于CFD网格的逐渐稀疏, 会导致一部分诱导阻力转化为熵增阻力。这就会使得诱导阻力偏小, 黏性阻力偏大, 但总阻力保持不变。因此, 远场法中积分区域的大小明显影响诱导阻力的大小。在具体计算过程中, 需要研究积分区域大小对阻力的影响并考虑诱导阻力修正项, 从而削弱积分区域大小对诱导阻力的影响。对于诱导阻力修正项, 由于诱导阻力积分区域可选择整个流场, 但是远场的翼尖涡耗散会导致诱导阻力转化为熵增阻力, 从而使得诱导阻力计算结果偏小。所以, 本文在计算诱导阻力时将涡区域内所产生的熵增阻力作为修正项添加到诱导阻力中, 以消除远场粗网格的数值耗散影响。

本文对于诱导阻力的积分面选择为距物体壁面一定距离以内的网格单元, 从图 9中可以看出诱导阻力数值随积分面选取不同的具体变化, 随着积分面距离的增加, 诱导阻力系数的预测结果变化较小, 当取距物体壁面一倍平均气动弦长的距离来进行积分就可取得较好的结果。因此最终决定在本文中将远场积分诱导阻力的区域选为, 距物体壁面一倍平均气动弦长以内的网格单元。

3 流场求解方法及求解器验证

本文采用基于三维可压缩RANS方程的求解器进行数值模拟求解, 其守恒形式的表达式为

(23)

式中:U为守恒量；F_{i(i=1, 23)}为无黏通量; G_{i(i=1, 2, 3)}为黏性通量。空间离散化方法为有限体积法, 时间推进方式为Directional Alternative Diagonal Implicit方法进行半隐式时间推进和Newton-Krylov方法进行全隐式时间推进^[18], 湍流模型采用S-A模型。

4 基于阻力分解方法的气动特性研究 4.1 网格收敛性研究

本小节将针对CRM翼身组合体构型开展网格收敛性研究。计算状态为Ma=0.5, C_L=0.5, Re=5×10⁶。所采用的1族5套计算网格均来自于DPW组委会官方, 具体的网格量和Y⁺值如表 1所示。

基于近场法和远场法在不同网格量下的气动力计算结果分别如表 2和表 3所示。

图 10给出了CRM构型基于近场法和远场法在不同网格量下的总阻力预测对比, 以及各项阻力分量占比的对比。此外, 图 11给出了构型在L3计算网格下的激波阻力区域和黏性阻力区域示意图。因此可以得到如下结论:

1) 从表 2可以看出, 基于近场法预测阻力, 计算网格从L1到L5, 阻力减小量是14.6 counts, 且从图 10b)可看出其中压差阻力减小量是13.1 counts, 而摩擦阻力的减小量只有1.5 counts。

2) 从表 3可以看出, 基于远场法预测阻力, 计算网格从L1到L5, 阻力减小量是10.4 counts, 相比近场法阻力的变化量减小了4.2 counts, 两者的具体对比如图 10a)所示; 而从图 10c)中可以看出远场法下各阻力分量的预测对比。

3) 从L1到L5, 黏性阻力减小量是21.7 counts, 其中黏性压差阻力减小量是20.2 counts, 比压差阻力的减小量大7.1 counts。

4) 从L1到L5, 可以发现激波阻力对总阻力的贡献只是很小一部分, 而且其递增变化的幅度很小, 只有1.2 counts。此外, 如果不考虑L1网格的话, 其变化的幅度更小, 只有0.3 counts。

5) 从L1到L5, 诱导阻力变化趋势和激波阻力情况类似, 也呈现增长趋势, 递增量为10.1 counts, 且从L4到L5的增长幅度只有1.4 counts。

6) 从伪阻力来看, 最大值出现在L1网格中, 为4.4 counts, 而最小值仅有0.2 counts, 是在最好的网格L5中得到的。而伪阻力从L1到L5的减小趋势也反映出网格质量的提升。此外, 伪阻力4.2 counts的递减量占压差阻力递减量的32%。

7) 从基于L5网格的计算数据可以看出:黏性阻力占总阻力的65%(而黏性阻力中70%来源于摩擦阻力, 30%为黏性压差阻力), 激波阻力仅占总阻力的2%, 而诱导阻力占全机总阻力的33%。

4.2 马赫数和雷诺数对气动力影响的分析研究 4.2.1 马赫数对气动力影响

本小节将针对CRM构型开展马赫数对气动力影响的分析研究。计算网格选取4.1节的L3网格, 所采用的计算状态为:C_L=0.5, Re=5×10⁶, Ma=0.70, 0.75, 0.80, 0.83, 0.85, 0.86, 0.87。

基于近场法和远场法在不同马赫数下的气动力计算结果如表 4所示。图 12为CRM构型基于近场法和远场法在不同马赫数下的总阻力预测对比, 以及各项阻力分量占比的对比。此外, 图 13给出了构型在不同马赫数下的激波阻力区域示意图。从这些数据对比可以总结如下结论:

1) 基于近场法预测阻力, 马赫数0.70~0.87, 阻力增加量是36.8 counts; 且从图 12b)中可以看出, 其中压差阻力从马赫数0.70时的124.2 counts增加到了马赫数0.85时的168.2 counts, 摩擦阻力呈递减趋势且变化量为7.2 counts, 这可能与机翼表面的流动分离有关。

2) 基于远场法预测阻力, 马赫数从0.70~0.87, 阻力递增量为36.2 counts, 相比近场法阻力的变化量减小了0.6 counts; 图 12c)给出了远场法下各阻力分量的对比, 其中激波阻力增加了21.8 counts, 约占压差阻力变化量的49.5%, 黏性压差阻力增加了20.6 counts, 约占压差阻力变化量的46.8%, 而诱导阻力和伪阻力分别增加了1 counts和0.6 counts。

3) 激波阻力从Ma=0.70时占远场法总阻力的1.6%左右, 增加到了Ma=0.85时占远场法总阻力的8.0%左右, 这从图 13激波区域的变化也可看出, 但激波阻力在总阻力中依然占比较小。

4.2.2 雷诺数对气动力影响

本小节将针对CRM翼身组合体构型开展雷诺数对气动力影响的分析研究。计算网格选取4.1节的L3网格, 所采用的计算状态为:C_L=0.5, Ma=0.85, Re=5×10⁶, Re=20×10⁶, Re=40×10⁶。

基于近场法和远场法在不同雷诺数下的气动力计算结果如表 5所示。图 14为CRM翼身组合体构型基于近场法和远场法在不同雷诺数下的总阻力预测对比, 以及各项阻力分量占比的对比。从这些数据对比可以总结如下结论:

1) 基于近场法预测阻力, 雷诺数从5×10⁶~40×10⁷, 阻力的减小量是37 counts; 且从图 14b)中可以看出其中摩擦阻力从Re=5×10⁶时的114 counts减小到了Re=40×10⁶时的90.6 counts, 减幅为23.4 counts, 约占近场法总阻力变化量的63%, 也说明了雷诺数的增加主要带来的就是边界层的变化进而引起摩擦阻力的变化, 而压差阻力也呈递减趋势且变化量为13.6 counts。

2) 基于远场法预测阻力, 雷诺数从Re=5×10⁶到Re=40×10⁶, 阻力的减小量是38.1 counts; 相比近场法阻力的变化量增大了1.1 counts, 这是由于伪阻力增加而造成的, 两者的具体对比如图 14a)所示; 图 14c)给出了远场法下各阻力分量的对比, 具体来说, 激波阻力减小了2.1 counts, 约占压差阻力变化量的15.4%, 黏性压差阻力减小了11.5 counts, 约占压差阻力变化量的84.6%, 说明了雷诺数的增加引起的压差阻力变化主要是因为边界层的改变进而引起尾迹区的变化, 而诱导阻力呈递减趋势, 但变化量较小, 为1.1 counts。

3) 综上来看, 雷诺数对阻力的影响主要体现在黏性压差阻力和摩擦阻力的改变, 二者变化量之和占到了远场法总阻力变化量的91.6%。

5 基于阻力分解的CRM构型气动外形优化设计

通常, 气动优化设计为了减小优化的计算代价, 多运用网格量较小的粗化网格。此时, 优化收益就会在一定程度上承担粗网格带来的伪阻力影响。而阻力分解方法既可以有效辨别出优化结果中物理机制下的减阻效果, 也可以分辨出伪阻力的占比变化, 对优化减阻的分析和实际达到的减阻有效性和可靠性有很大帮助。

5.1 优化设计系统

本文的气动外形优化设计系统^[19-20]主要由外形参数化、网格变形算法、流场求解以及优化算法等模块组成, 已经由作者和课题组其他学者应用于相关民用客机的气动优化设计当中^[19-20]。

5.2 优化问题介绍

本文选取CRM翼身组合体构型作为气动优化设计研究的对象, 其FFD控制体如图 15所示。

从图中可以看出, 控制机翼作为设计变量的FFD控制点共有420个(蓝色点), 其中沿展向分布10个, 沿弦向分布21个, 垂直于弦向2个。而优化设计过程中采用的计算网格为4.1节中的L3网格, 网格量为5 111 808, 基于此本文开展了涉及大规模设计变量考虑多约束的气动减阻优化设计。优化的巡航设计点状态是C_L=0.5, Ma=0.85, Re=5×10⁶, 且考虑到工程实际应用, 跨声速民用宽体客机气动设计的重要指标之一是具有良好的阻力发散特性, 因此, 优化的另一个设计工况为C_L=0.5, Ma=0.87, Re=5×10⁶。优化目标是巡航设计点的阻力系数和阻力发散马赫数下的阻力系数加权后总阻力系数最小, 优化设计变量是420个机翼型面的FFD控制点和2个工况下的自由来流攻角(实现定升力系数)。此外, 由于机翼需满足一定结构要求和容积要求, 因此必须对机翼的厚度进行约束。本文在优化中采取了150个厚度约束, 其中沿展向分布10个, 弦向均布15个, 保证每一个约束处的相对厚度t_i不小于初始值t_{i_initial}。综上, 该优化问题的数学模型即为

(24)

5.3 优化设计结果及分析

CRM翼身组合体构型优化前后的气动特性对比如表 6所示(其中Ma0.85_ini和Ma0.87_ini分别表示初始构型在Ma=0.85下和Ma=0.87下的计算结果, Ma0.85_opt和Ma0.87_opt分别表示优化构型在Ma=0.85下和Ma=0.87下的计算结果)。图 16为构型优化前后在近场法和远场法下的总阻力预测对比, 以及各项阻力分量占比的对比。图 17给出了构型优化前后在不同马赫数下的激波阻力区域示意图。从这些数据对比可以看出:

1) 从表 6中基于近场法的结果可以看出, 优化构型巡航设计点的总阻力减小了3.8 counts, 减幅为1.5%, 且升阻比从19.8提升到20.1;而阻力发散马赫数下的总阻力减小了23.1 counts, 减幅为8.2%, 且升阻比从17.8提升到19.4。

2) 从图 16b)中可以看出近场法下的阻力分量对比, 其中优化构型巡航设计点的压差阻力从138.1 counts减小到了134.2 counts, 摩擦阻力增加了0.1 counts; 阻力发散马赫数下的压差阻力从168.2 counts减小到了144.2 counts, 减幅非常明显, 摩擦阻力增加了0.9 counts; 2个工况下的摩阻增加均可能与机翼表面激波强度的改变有关。

3) 构型经过优化设计后, 从表 6中基于远场法预测的力系数可以看出, 巡航设计点的总阻力减小了3.7 counts, 较近场法减阻量少了0.1 counts, 但总阻力较近场法减少了0.7 counts, 两者的具体对比如图 16a)所示, 这也体现出来利用远场法去除伪阻力之后流场中基于物理机制成因达到的减阻量, 且升阻比从19.9提升到20.2, 优化后升阻比较近场法下又增加了0.1;而阻力发散马赫数下的总阻力减小了21.9 counts, 虽然较近场法减阻量减少了1.2 counts, 但总阻力较近场法下减少了0.5 counts, 另外升阻比从17.9提升到19.5, 优化后升阻比较近场法下增加了0.1。

4) 图 16c)给出了远场法下的各阻力分量对比, 这也为近场法下优化前后压差阻力的减小提供了更加详细的阻力占比分析。机翼经过优化设计后, 巡航设计点的激波阻力减小了3.8 counts, 黏性压差阻力增大了1.0 counts, 且诱导阻力和伪阻力分别减小了1.0 counts和0.1 counts; 而阻力发散马赫数下的激波阻力减小了15.7 counts, 占到此工况下减阻量的71.7%, 黏性压差阻力减小了5.5 counts, 占到此工况下减阻量的25.1%, 也是减阻的一部分重要贡献, 且诱导阻力和伪阻力分别减小了1.6 counts和1.2 counts。

综上可以说明, 将改进的阻力分解方法融入基于梯度的气动优化设计系统, 对于三维构型考虑多工况的气动优化设计, 不仅可以更加有效分析减阻的来源和对应阻力区域的变化, 还可以进行不同工况下的阻力变化对比, 对于优化前后构型的气动性能变化有更加深刻的理解, 这也为设计师提供了更加直观可靠的数据和性能分析。

6 结论

1) 为了能够准确可靠的进行阻力预测且有效提升阻力预测的精度, 针对阻力预测方法开展了相关研究。首先针对国内外现有的几种利用远场法进行阻力预测和分解的方法进行了对比讨论分析, 推导了相关理论公式和发展方向, 提炼出各自方法针对轴向速度损失量的优劣势和差异。进而建立了本文关于轴向速度损失量的改进方法, 提出了本文基于远场法的阻力预测方法和阻力分解方法。

2) 为了能够更加准确可靠地使用本文所推导建立的阻力分解方法, 针对阻力区域的判断选择进行了详细的参数敏感性分析研究。首先, 构造了一种区域划分算法, 判定不同阻力区域识别的前后顺序; 其次, 针对不同的阻力进行区域判定, 选择了合适的激波区域判断识别函数和边界层尾迹区区域判断识别函数; 然后, 针对识别函数不能完全涵盖阻力选择区域的问题, 进行了阻力区域的扩增判断和识别函数值的阈值选定; 此外, 基于诱导阻力的积分定义, 对其积分面的选取也进行了分析和确定。

3) 基于阻力分解方法, 针对CRM翼身组合体构型, 进行了详细的气动特性研究。首先进行了网格收敛性研究, 结果显示随着网格从L1到L5的改善和提升, 基于近场法的总阻力系数递减量为14.6 counts, 而基于远场法的递减量可以减小4.2 counts, 且不同阻力分量大小和其区域可以进行直观的对比, 以及力系数的预测精度也可以得到充分保证; 其次, 马赫数和雷诺数的影响研究展现出更丰富详细的阻力变化对比, 且通过不同阻力分量的变化对比有效的阐述和解释了其对总阻力的变化影响以及具体的贡献占比。

4) 将改进的阻力分解方法融入基于梯度的气动外形优化设计系统, 针对CRM翼身组合体构型进行了气动外形优化设计研究。优化设计兼顾了巡航设计点和阻力发散工况点, 且涉及大规模设计变量和约束条件。优化结果详细对比分析了三维复杂构型在多工况条件下的减阻量对比以及不同阻力分量变化对于减阻的占比贡献, 而且通过阻力区域识别函数直观感受阻力分量变化的详细情况, 此外, 可以更加精确的得到优化构型去除伪阻力以后的总阻力与升阻比。