对于在超音速飞行状态下的飞行器，激波阻力占据总阻力的很大一部分。1935年，Adolf Busemann[1]提出了Busemann双翼的概念。将2个三角翼反对称放置，利用激波和膨胀波的相互作用，在理论上Busemann双翼可以完全消除波阻。Kusunose等[2]使用计算流体力学(CFD)技术验证了Busemann双翼具有良好的减阻特性。华如豪等[3]将Busemann双翼应用于超音速导弹，发现在巡航条件下，双翼布局能够使波阻降低42%。然而，Busemann双翼是在特定条件下设计的。当双翼偏离设计工作状态时，气动特性有可能急剧变化，会存在流动壅塞问题，导致阻力急剧升高[2]。因此，有必要分析来流参数不确定性对双翼气动特性的影响，寻找出双翼气动特性波动的主要原因，为双翼的稳健性优化设计提供参考。
传统的气动布局设计仅仅考虑在确定状态下的气动特性，而往往忽略不确定性参数带来的影响[4-5]。近年来，随着CFD技术的迅速发展，人们更加追求计算结果的精度及可信度[6]。因此，不确定性分析方法的研究及应用受到了广泛关注。蒙特卡罗方法由于结构简单，易于实现，较早地应用在不确定性分析方法中。但是由于需要大量样本，计算开销大，蒙特卡罗方法在实际应用中存在一定困难。近年来，基于多项式混沌(PC)的不确定性分析方法则逐渐发展起来[7-8]。多项式混沌可以分为嵌入式多项式混沌(IPC)和非嵌入式多项式混沌(NIPC)。IPC方法比较直观，对控制方程进行混沌多项式展开，进而可以得到相关的统计信息。但是由于改变了控制方程，需要对求解器进行大量更改，工作量较大[9]。NIPC方法是将通过抽样方法得到的样本点直接带入求解器中，从而获取确定性的解。依据这些确定性的结果，通过回归分析的方法，求解混沌多项式的各项系数，得到近似的混沌多项表达式。相比于IPC方法，NIPC方法不需要对求解器进行修改，降低了工作量。但是在高维问题中，依然存在工作量较大的问题[10]。在保证结果精度的前提下，为了进一步减少计算时间，Loeven等[11]提出了非嵌入式概率配置点(NIPRC)方法。NIPRC方法是在NIPC方法的基础上，使用配置点代替抽样样本点。Loeven[12]将IPC, NIPC及NIPRC方法应用于一维线性活塞问题。结果表明，NIPRC方法比NIPC方法收敛更快。在相同精度时，NIPRC方法比IPC方法计算耗时更小。Dinescud等[13]验证了NIPRC和NIPC方法在NASA rotor转子算例中结果具有较好的一致性。刘益智等[14]使用NIPRC方法分析了转子间隙不确定性对风力机气动特性的影响，发现前缘间隙的不确定性是引起气动性能变化的决定性因素。张航等[15]应用NIPRC方法研究后台阶流动。
单纯的不确定性分析只能获得不确定性因素对输出的影响，为了进一步筛选出影响显著的变量，这就需要对模型进行敏感性分析[16-17]。敏感性分析分为局部敏感性和全局敏感性[18]。局部敏感性可以计算单一变量对模型的影响程度。全局敏感性不仅可以衡量各个变量对模型的贡献，也可以求得各个变量之间耦合的影响。常用的全局敏感性分析方法有多元线性回归、傅里叶幅度敏感性检验和Sobol等方法。其中，基于方差的Sobol方法[19]具有计算迅速和可操作性强等优点，是最具代表性的全局敏感性方法。
本文使用NIPRC方法对二维Busemann双翼气动特性进行不确定性量化研究。分析了在来流马赫数、迎角等不确定因素作用下，双翼流场的变化情况，得到造成双翼气动特性波动的主要原因。并进一步使用Sobol全局敏感性方法，寻找出对双翼气动特性影响显著的因素，分析马赫数和迎角耦合作用对双翼升力和阻力的影响。
1 非嵌入式概率配置点方法依据非嵌入式概率配置点思想，输出函数是由确定和随机两部分组成的。对于任意一个属于样本空间Ω上的随机变量u(x, t, ω), 可以用下式逼近表示[12]
|
(1) |
式中, u(x, t, ω)为空间x和时间t的随机变量, ω为样本空间Ω的基本事件, Np为配置点的个数, 系数ui(x, t)为依据配置点得到的确定性部分。Li是与配置点相对应的拉格朗日多项式混沌, 其表达式如下
|
(2) |
对于数值型求积问题, 如果n+1个求积节点给定, 仅选取求积系数, 使用Gauss型求积公式代数精度最高, 为2n+1。将选取的高斯积分点映射到样本空间, 可得到配置点。对于服从正态分布的随机变量, 配置点可以通过Gauss-Hermite求积公式得到。通过查表得到Gauss-Hermite求积节点和系数, 并将其映射到样本空间, 即可得到配置点及相应权重。
通过(1)式, 不难获得输出函数的相关统计信息。输出函数的平均值、方差和变异系数可以表示为
|
(3) |
|
(4) |
|
(5) |
式中, ωi是配置点相应的权重。
由于非嵌入式配置点方法只选取了有限阶模态展开, 必然存在误差。相对误差可由更高一阶的表达式得到[21]。P阶相对误差的表达式为
|
(6) |
式中, ωi是P+1阶配置点对应权重。uiP+1是P+1阶配置点对应的确定性部分。
全局灵敏度分析不仅可以检验模型各个参数变化对输出的影响, 还可以分析各个参数之间耦合作用对输出的影响。Sobol是基于方差的全局敏感性分析方法[19]。总方差可以表示为
|
(7) |
式中，Vi=V[E(Y|Xi)], Vij=V[E(Y|Xi, Xj)], …, 将上式两边除以V(Y), 得到
|
(8) |
式中, Si为一阶灵敏度指标, Sij…k为高阶灵敏度指标。Sij…k定义为
|
(9) |
变量Xi总灵敏度指标为
|
(10) |
本文采用Saltelli[22]发展的Sobol全局灵敏度计算方法:
1) 产生一个(N, 2K)随机变量矩阵, 其中K为输入变量的个数, N为抽取样本数。将矩阵分为两个矩阵A和B
|
2) 在B矩阵的基础上, 将A矩阵的第i列与B矩阵的第i列交换, 得到Ci矩阵
|
3) 根据输入矩阵A, B, Ci计算相应的输出
|
4) 主灵敏度指标计算如下
|
(11) |
式中
|
(12) |
5) 总灵敏度指标计算如下
|
(13) |
6) 对于二维问题, 二阶灵敏度指标计算如下
|
(14) |
图 1展示了Busemann双翼的消波原理。其中c为弦长, t为单翼厚度, z为双翼间距。通过合理的参数设计, 使得2个三角形单元前缘产生的激波恰好打在翼型最大厚度处, 最大厚度处的膨胀波又恰好打到三角形后缘[3]。这样气流在激波和膨胀波作用下, 理论上在某一马赫数Busemann双翼迎风面和背风面压力几乎相等, 波阻接近为零。此时, Busemann双翼构型达到设计状态, 对应的马赫数称为设计马赫数。
|
| 图 1 Busemann双翼模型示意图 |
本文考虑的输入变量有2个, 分别是来流马赫数和迎角。变量之间相互独立, 使用NIPRC方法进行不确定性量化研究, 研究来流状态变量不确定性对双翼气动特性的影响。
4.1 流场求解及网格无关性检验本文采用课题组自研制基于有限体积法的流场求解程序HUNS3D进行流场计算[23]。HUNS3D求解程序的计算精度已得到验证[23-24]。空间离散格式为Roe格式, 时间推进格式为LU-SGS。选取NS方程作为控制方程。自由来流马赫数为1.7, 来流迎角为2°, 弦长c, 单翼厚度t=0.05c。双翼间距z=0.5c。如图 2所示, 采用结构网格划分计算区域。
|
| 图 2 计算网格 |
为确保计算所用网格的数量与结果之间没有关联性, 进行网格无关性检验。结果见表 1。可以看到, 当网格单元数量达到108 375时, 网格数量对双翼的升力、阻力系数的影响几乎可以忽略。因此, 选取该网格为计算网格。双翼中上单翼的压力分布如图 3所示。激波和膨胀波先后作用在上单翼最大厚度处, 这导致该区域的压强先陡然提升后又陡然下降。受到来流迎角的影响, 由下单翼反射出的激波作用于上单翼下表面距离前缘0.8倍弦长附近, 紧接着在该区域又反射出膨胀波。其具体表现为, 压强在该区域先提升后下降。
| 网格单元数量 | CL | CD |
| 46 151 | 0.117 | 0.0163 3 |
| 73 317 | 0.118 | 0.0164 1 |
| 108 375 | 0.119 | 0.0165 4 |
| 210 477 | 0.119 | 0.0165 7 |
|
| 图 3 上单翼的压力分布 |
假定马赫数为服从均值1.7, 标准差0.011 56正态分布的随机变量, 迎角为服从均值2, 标准差0.071 6正态分布的随机变量。马赫数和迎角的概率密度函数如图 4所示。使用NIPRC方法研究各个变量独立变化对双翼气动特性的影响。为了验证NIPRC方法的收敛性, 分别选取P=0, 1, 2阶配置点进行确定性计算, 按公式(6)计算相对误差。
|
| 图 4 马赫数和迎角的概率密度函数 |
| 多项式 阶数 |
配置点 | 权重 | 马赫数 | 迎角 |
| 0 | 1 | 1 | 1.7 | 2 |
| 1 | 1 2 |
0.5 0.5 |
1.688 44 1.711 56 |
1.928 40 2.071 60 |
| 2 | 1 2 3 |
0.166 67 0.666 67 0.166 67 |
1.679 98 1.700 00 1.720 02 |
1.875 99 2.000 00 2.124 01 |
| 3 |
1 2 3 4 |
0.045 88 0.454 12 0.454 12 0.045 88 |
1.673 01 1.691 42 1.708 58 1.726 99 |
1.832 86 1.946 88 2.053 12 2.167 14 |
图 5展示了随着NIPRC方法中多项式阶数的增加, 双翼升力、阻力系数的相对误差。当多项式阶数P=2时, 升力系数和阻力系数的相对误差在10-3量级。此时, 升力的绝对误差低于一个升力单位(一个升力单位为10-2), 阻力的绝对误差低于一个阻力单位(一个阻力单位为10-4)。因此, 本节选取二阶NIPRC方法做进一步分析。
|
| 图 5 马赫数不确定性条件下, 双翼气动特性随多项式阶数的收敛过程 |
使用二阶方法得到上单翼压强的平均值和标准差, 见图 6。标准差越大, 马赫数的随机扰动影响越强烈, 上单翼上表面压强的标准差接近0, 这表明马赫数的随机扰动对上单翼上表面压强影响极小。上单翼压强标准差最大的地方是下表面的中间区域(翼型最大厚度处)和下表面距离前缘0.7倍弦长附近。究其原因, 受到马赫数的不确定性影响, 上单翼下表面中间区域激波和膨胀波变化剧烈。而由下单翼反射的激波打到了距离上单翼下表面前缘0.7倍弦长附近。这导致该区域的压强有小幅波动。
|
| 图 6 马赫数不确定性条件下, 上单翼压强平均值和标准差 |
为进一步研究马赫数的不确定性对整个流场的影响, 图 7给出了全流场压强平均值和标准差云图。
|
| 图 7 马赫数不确定性条件下, 全流场压强平均值和标准差 |
由压强的平均值云图可以看出, 在马赫数随机扰动的影响下, 压强在下单翼上表面区域有明显的增加。而标准差数值较大的区域主要集中在激波和膨胀波处。其中, 下单翼上表面中间区域的压强标准差数值最大, 这说明该区域受到的影响程度最大。其表现为该区域压强产生较大波动, 从而严重影响双翼的气动特性。由图 8可知, 升力系数接近正态分布, 而阻力系数接近偏斜分布。表 3给出了在马赫数随机扰动下, 双翼的升、阻力系数的平均值和标准差。当输入参数马赫数的变异系数为0.68%时, 升力系数的变异系数为2.46%, 阻力系数的变异系数为2.61%。这说明马赫数的不确定性可引起双翼气动特性的强烈变化。
|
| 图 8 马赫数不确定性条件下, 双翼升力和阻力系数的概率密度函数 |
| 不确定参数 | CL | CD | ||||
| μ | σ | CV/% | μ | σ | CV/% | |
| 马赫数 | 0.119 3 | 0.002 93 | 2.46 | 0.016 61 | 0.004 33 | 2.61 |
| 迎角 | 0.119 0 | 0.004 35 | 3.66 | 0.016 55 | 0.000 32 | 1.93 |
| 确定性条件 | 0.118 8 | 0.016 54 | ||||
双翼升力、阻力系数的相对误差随着NIPRC方法阶数增加的收敛过程见图 9。当多项式阶数P=2时, 升力系数的相对误差在10-4量级, 阻力系数的误差在10-5量级。此时, 升力的绝对误差低于一个升力单位, 阻力的绝对误差低于一个阻力单位。本节选取二阶NIPRC方法做进一步分析, 二阶NIPRC方法需要3个配置点做确定性计算。
|
| 图 9 迎角不确定性条件下, 双翼气动特性随多项式阶数的收敛过程 |
上单翼压强的平均值和标准差图如图 10所示。从压强标准差图可知, 压强标准差的最大值为0.031。与之对比的是, 在马赫数扰动下, 上单翼压强标准差的最大值为0.082。标准差越小则受到的影响越小。所以, 相比于马赫数, 由迎角带来的不确定性对上单翼压强的波动贡献更小。
|
| 图 10 迎角不确定性条件下, 上单翼压强平均值和标准差 |
整个流场压强平均值和标准差云图见图 11。
|
| 图 11 迎角不确定性条件下, 全流场压强平均值和标准差 |
结果表明, 流动变化剧烈的区域主要集中在下单翼上表面翼型最大厚度处。由图 8可知, 升力和阻力系数的概率分布趋势接近高斯分布。表 3给出了在迎角随机扰动下, 双翼的升、阻力系数的平均值与标准差。当输入参数迎角的变异系数为3.58%时, 升力系数的变异系数是3.66%, 阻力系数的变异系数为1.93%。这说明双翼的气动参数对迎角的变化并不敏感。
|
| 图 12 迎角不确定性条件下, 双翼升力和阻力系数的概率密度函数 |
表 3给出了在马赫数和迎角各自独立随机扰动下, 双翼升力、阻力系数不确定性结果。从输出结果不难发现, 马赫数比迎角更易引起双翼的气动特性波动。当不确定参数为迎角时, 双翼升力系数的变异系数高于阻力系数的变异系数。这说明升力比阻力更易受到迎角扰动的影响。而马赫数对阻力的影响高于对升力的影响。为了进一步衡量不确定参数马赫数和迎角的影响程度并分析马赫数和迎角耦合作用对双翼气动特性的影响, 需要进行Sobol全局敏感性分析。首先通过拉丁超立方方法随机抽取1 000个样本点。其次, 调用流场数值模拟工具HUNS3D进行确定性计算, 得到各样本点对应的输出结果(升力和阻力系数)。最后, 依据第2节所述的计算方法, 得到全局灵敏度指标, 结果见表 4。其中, SM为马赫数的一阶灵敏度指标。Sα为迎角的一阶灵敏度指标。SM-α为二阶灵敏度指标, TM代表了马赫数总灵敏度, Tα代表了迎角总灵敏度。结果表明:①升力系数中迎角的一阶灵敏度指标最高, 这说明迎角对双翼的升力变化占主要贡献。阻力系数中马赫数的一阶灵敏度指标最高, 这说明马赫数对双翼的阻力变化占主要贡献。②变异系数为0.68%的马赫数对升力和阻力系数的总灵敏度指标分为0.459和0.591, 变异系数为3.58%的迎角对升力和阻力系数的总灵敏度指标分别为0.615和0.416, 这说明造成双翼气动特性波动最主要的因素是马赫数; ③马赫数和迎角对升力和阻力的二阶灵敏度指标为0.074和0.043, 这说明马赫数和迎角的耦合作用对升力和阻力波动的贡献较小, 可以忽略。
| 灵敏度指标 | CL | CD |
| SM | 0.385 | 0.584 |
| Sα | 0.541 | 0.373 |
| SM-α | 0.074 | 0.043 |
| TM | 0.459 | 0.591 |
| Tα | 0.615 | 0.416 |
本文采用非嵌入式概率配置点不确定性分析方法和Sobol全局灵敏度分析方法，分析了马赫数和迎角的不确定性对双翼气动特性的影响，结论如下：
1) 不确定性分析表明，迎角对双翼升力的影响高于对阻力的影响。马赫数不确定性的影响区域主要集中在双翼间激波和膨胀波处。其中，下单翼上表面翼型厚度最大处受到的影响最严重。这导致该区域的压强脉动，从而严重影响双翼的气动性能。相比于马赫数，迎角的不确定性影响区域更小。
2) 灵敏度分析定量给出了迎角和马赫数的不确定性对双翼气动特性的影响程度。其中，马赫数是造成双翼气动特性波动的最主要因素。而马赫数和迎角的耦合作用对升力和阻力的影响相对较小。
| [1] | BUSEMANN A. Aerodynamic Lift at Supersonic Speeds[J]. Luftfahrtforschung, 1935, 12(6): 210-220. |
| [2] | KUSUNOSE K, MATSUSHIMA K, MARUYAMA D. Supersonic Biplane-A Review[J]. Progress in Aerospace Sciences, 2011, 47(1): 53-87. DOI:10.1016/j.paerosci.2010.09.003 |
| [3] |
华如豪, 叶正寅. 基于Busemann双翼构型的超音速导弹减阻技术研究[J]. 应用力学学报, 2012, 29(5): 535-540.
HUA Ruhao, YE Zhengyin. Drag Reduction Method for Supersonic Missile Based on Busemann Biplane Concept[J]. Chinese Journal of Applied Mechanics, 2012, 29(5): 535-540. (in Chinese) |
| [4] | LI W, HUYSE L, PADULA S. Robust Airfoil Optimization to Achieve Drag Reduction over a Range of Mach Numbers[J]. Structural & Multidisciplinary Optimization, 2002, 24(1): 38-50. |
| [5] |
高正红, 王超. 飞行器气动外形设计方法研究与进展[J]. 空气动力学学报, 2017, 35(4): 516-528.
GAO Zhenghong, WANG Cao. Aerodynamic Shape Design Methods for Aircraft:Status and Trends[J]. Acta Aerodynamica Sinica, 2017, 35(4): 516-528. (in Chinese) DOI:10.7638/kqdlxxb-2017.0058 |
| [6] | SLOTNICK J, KHODADOUST A, ALONSO J, et al. CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences[R]. NASA/CR-2014-218178 |
| [7] | DOOSTAN A, GHANEM R G, RED-HORSE J. Stochastic Model Reduction for Chaos Representations[J]. Computer Methods in Applied Mechanics & Engineering, 2007, 196(37/38/39/40): 3951-3966. |
| [8] | MATRE O P L, KNIO O M. Spectral Methods for Uncertainty Quantification[M]. Berlin: Springer Netherlands, 2010. |
| [9] | XIU D, KARNIADAKIS G E. Modeling Uncertainty in Flow Simulations Via Generalized Polynomial Chaos[J]. Journal of Computational Physics, 2003, 187(1): 137-167. |
| [10] | HOSDER S, WALTERS R. A Non-Intrusive Polynomial Chaos Methods for Uncertainty Propagation in CFD Simulation[C]//44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nedava, 2006: 1-19 |
| [11] | LOEVEN G J A, BIJL H. Probabilistic Collocation Used in a Two-Step Approached for Efficient Uncertainty Quantification in Computational Fluid Dynamics[J]. Computer Modeling in Engineering & Sciences, 2008, 36(3): 193-212. |
| [12] | LOEVEN G J A. Efficient Uncertainty Quantification in Computational Fluid Dynamics[D]. Holland, Delft University of Technology, 2010 |
| [13] | DINESCU C, SMIRNOV S, HIRSCH C. Assessment of Intrusive and Non-Intrusive Non-Deterministic CFD Methodologies Based on Polynomial Chaos Expansions[J]. International Journal of Engineering Systems Modelling & Simulation, 2010, 2(1): 87-98. |
| [14] |
刘智益, 王晓东, 康顺. 叶顶间隙尺度的不确定性对压气机性能影响的CFD模拟[J]. 工程热物理学报, 2013, V34(4): 628-631.
LIU Zhiyi, WANG Xiaodong, KANG Shun. CFD Simulations of Uncertain Tip Clearance Effect on Compressor Performace[J]. Journal of Engineering Thermophysics, 2013, V34(4): 628-631. (in Chinese) |
| [15] |
张航, 吴新. 非嵌入式概率配置点法在后台阶地形流场计算中的应用[J]. 水电能源科学, 2017(7): 30-32.
ZHANG Han, WU Xin. Application of Non-Intrusive Probabilistic Collocation Method in Back-Step Flow Calcultion[J]. Water Resources and Power, 2017(7): 30-32. (in Chinese) |
| [16] |
宋赋强, 阎超, 马宝峰, 等. 锥导乘波体构型的气动特性不确定度分析[J]. 航空学报, 2018, 39(2): 121519.
SONG Fuqiang, YAN Chao, MA Baofeng, et al. Uncertainty Analysis of Aerodynamic Characteristics for Cone-Derived Waverider Configure[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(2): 121519. (in Chinese) |
| [17] |
邬晓敬, 张伟伟, 宋述芳, 等. 翼型跨声速气动特性的不确定性及全局灵敏度分析[J]. 力学学报, 2015, 47(4): 587-595.
WU Xiaojing, ZHANG Weiwei, SONG Shufang, et al. Uncertainty Quantification and Global Sensitivity Analysis of Transonic Aerodynamics about Airfoil[J]. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(4): 587-595. (in Chinese) |
| [18] |
蔡毅, 邢岩, 胡丹. 敏感性分析综述[J]. 北京师范大学学报, 2008, 44(1): 9-16.
CAI Yi, XING Yan, HU Dan. On Sensitivity Analysis[J]. Journal of Beijing Normal University, 2008, 44(1): 9-16. (in Chinese) DOI:10.3321/j.issn:0476-0301.2008.01.003 |
| [19] | SOBOL I M. Sensitivity Estimates for Nonlinear Mathematical Models[J]. Math Model Comput Exp, 1993, 1(1): 112-118. |
| [20] | SOBOL I M. Global Sensitivity Indices for Nonlinear Mathematical Models and Their Monte Carlo Estimates[J]. Mathematics & Computers in Simulation, 2014, 55(1): 271-280. |
| [21] | TATANG M A, PAN W, PRINN R G, et al. An Efficient Method for Parametric Uncertainty Analysis of Numerical Geophysical Models[J]. Journal of Geophysical Research Atmospheres, 1997, 102(18): 21925-21932. |
| [22] | Saltelli. Global Sensitivity Analysis[M]. England: John Wiley, 2008. |
| [23] |
王刚, 叶正寅. 三维非结构混合网格生成与N-S方程求解[J]. 航空学报, 2003, 24(5): 385-390.
WANG Gang, YE Zhengyin. Generation of Three Dimensional Mixed and Unstructured Grids and its Application in Solving NavierStokes Equations[J]. Acta Aeronautica et Astronautica Sinica, 2003, 24(5): 385-390. (in Chinese) DOI:10.3321/j.issn:1000-6893.2003.05.001 |
| [24] | WANG G, YE Z Y, WANG G, et al. Application of Mixed Element Type Unstructured Grid in Solving Navier-Stokes Equations[J]. Chinese Journal of Computation Physics, 2004, 21(2): 161-165. |
