2. 中国空气动力研究与发展中心空气动力学国家重点实验室，四川绵阳 621000
2. State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000, China
2009年6月召开的第4届DPW会议(DPW IV)选择了CRM翼/身/平尾组合体构型(common research model wing-body-horizontal tail configuration, CRM-WBH)作为基准研究模型[1], 该模型分别在NASA Langley的NTF风洞(National Transonic Facility), NASA Ames的11英尺TWT风洞(Transonic Wind Tunnel)和欧洲的ETW风洞(European Transonic Wind Tunnel)开展了风洞试验, 试验结果包括了气动特性、表面压力分布及模型变形测量数据等[2-3], 这些公开发布的试验结果为CFD的验证和确认工作提供了丰富的基础数据. DPW IV多家单位采用CFD方法得到的气动特性计算结果与试验结果的对比分析显示[4], 数值模拟得到的相同迎角下的升力系数普遍大于试验结果, 俯仰力矩系数普遍低于试验结果, 尤其是俯仰力矩系数与试验结果差别明显.计算模型中没有考虑的试验模型的支撑装置和气动载荷作用下的静气动弹性变形是导致上述差异的主要原因.
2012年, Rivers等[5-6]利用NTF风洞试验测量得到的模型变形数据, 采用非结构网格技术和CFD方法研究了风洞试验模型的支撑机构和静气动弹性变形对CRM-WBH构型数值模拟结果的影响, 但此项工作只采用了升力系数等于0.5时的变形测量数据, 而没有考虑试验模型的静气动弹性变形量随来流迎角的变化. 2015年, Keye等[7]采用非结构网格技术和流固耦合方法模拟了CRM-WBH构型的气动特性, 但此项工作中没有考虑模型支撑装置对数值模拟结果的影响.
利用DPW组委会提供的CRM-WBH计算模型、结构有限元模型和支撑模型, 基于结构网格技术, 采用CFD方法和流固耦合计算方法开展了包含支撑装置的CRM-WBH模型(CRM-WBHS)数值模拟, 通过与ETW风洞测力、测压和模型变形测量结果的对比分析(试验数据取自http://common research model.larc.nasa.gov/experimental approach), 研究了模型支撑装置和静气动弹性变形对CRM-WBH模型气动特性数值模拟结果的影响.
1 CRM-WBH计算模型与气动计算网格CRM模型是典型的现代运输机构型, 巡航Mach数为0.85, 名义升力系数为0.50, 设计外形的详细参数见文献[1], 本文选择了CRM-WBH构型作为研究对象, 平尾的安装角为0°. CRM模型试验外形模型缩比为0.027, 基本参数如下:模型参考面积Sref = 0.279 7 m2, 平均气动弦长c = 0.189 m, 展长b = 1.587 m，梢根比λ = 0.275, 展弦比AR = 9.0, 1/4弦线后掠角Λc/4 = 35.0°.计算模型见图 1.
|
| 图 1 CRM-WBH计算模型 Fig.1 CRM-WBH simulation model |
CRM-WBH风洞试验模型采用安装于机身后体的叶片尾撑方式固定于风洞迎角变换装置.由于风洞试验模型后部的迎角变换装置对气动特性的影响很小[5], 因此本文的数值模拟中, 没有考虑模型支撑后部的迎角变换装置; 同时, 对模型叶片尾撑延伸段进行了局部修型处理以避免底部分离导致的计算收敛困难.本文将经过上述处理的带支撑装置的CRM-WBH模型标识为CRM-WBHS(见图 2).
|
| 图 2 CRM-WBHS计算模型(局部) Fig.2 CRM-WBHS simulation model(local) |
NASA Langley研究中心在互联网上公布了一套采用四面体实体单元离散的CRM风洞模型结构有限元模型, 整个有限元模型包含约1.4×106个网格节点, 6.8×106个网格单元和8.2×106个自由度.为与本文的计算构型和风洞试验模型相一致, 本文采用的CRM-WBH结构有限元模型去除了原始有限元模型中的挂架和短舱(见图 3).
|
| 图 3 CRM-WBH模型有限元模型 Fig.3 CRM-WBH finite element model |
本文采用CFD方法和流固耦合(fluid structure coupling, FSC)方法开展了CRM机翼/机身/平尾组合体气动特性的数值模拟, 以下简要介绍两种计算方法.
CFD计算方法:采用中国空气动力研究与发展中心研发的亚跨超CFD软件平台(TRIsonic Platform, TRIP), TRIP软件经过了系统的验证和确认工作[8-9].在本文的研究工作中, Reynolds平均Navier-Stokes方程无黏项离散采用2阶精度MUSCL格式[10], 黏性项离散采用2阶中心格式, 湍流模型采用Menter的SST两方程模型[11], 离散方程组的求解采用LU-SGS方法[12], 采用多重网格技术和大规模并行技术加速收敛.本文以下的数值计算中, 采用“全湍流”模拟方式, 没有考虑风洞试验的转捩位置.
流固耦合计算方法:流固耦合计算方法主要包含计算流体力学求解模块(CFD)、计算结构力学(computational structural mechanics, CSM)求解模块、耦合界面数据传递和动态网格变形4个主要功能模块, 采用松耦合方式有序组织上述功能模块, 实现复杂飞行器静气动弹性数值模拟.其中, CFD求解模块, 采用TRIP软件; CSM计算模块, 采用柔度矩阵方法求解线性结构静力学方程, 获得气动载荷作用下的物面变形; 流固耦合界面数据传递模块, 采用薄板样条(thin plate spline, TPS)插值方法[13]构建CFD计算模块与CSM结构计算模块之间的气动载荷与结构变形传递矩阵; 动态网格变形模块, 采用基于径向基函数(radial basis functions, RBF)[14]与超限插值(transfinite interpolation, TFI)[15]的复合型动态网格变形方法RBF_TFI.本文以下的流固耦合计算中, 只考虑了机翼变形, 机身和尾翼均采用刚体假设.
3 CRM-WBH模型网格收敛性研究根据DPW组织委员会提供的网格生成指导原则, 采用ICEM(Integrated Computer Engineering and Manufacturing Code)软件生成了CRM-WBH构型的极粗、粗、中、细不同网格规模的对接结构网格以开展网格收敛性研究, 为以下的数值模拟提供基础网格. 4套网格的详细信息参见表 1, 其中, Nnode表示网格节点总数, nBL, λBL分别表示边界层网格数量和网格增长率, y+为第1层物面网格法向无量纲距离. 图 4给出了CRM-WBH模型计算网格拓扑和表面网格(中等).
|
下载CSV 表 1 CRM-WBH构型网格参数 Tab.1 Grid parameters of CRM-WBH model |
|
| 图 4 CRM-WBH模型网格拓扑及表面网格(中等) Fig.4 Grid topology and surface grids of CRM-WBH model (medium grid) |
采用第2节介绍的CFD计算方法开展了固定迎角下的网格收敛性研究.计算来流条件为: Ma = 0.85, 迎角α = 3.0°, Re = 5.0×106. 表 2给出了采用不同密度网格计算得到的CRM-WBH模型的升力系数(CL), 阻力系数(CD)和俯仰力矩系数(Cm), 由表 2看出, CL, CD随网格密度的增加而单调变化, 从中等网格到密网格, CL变化0.000 2, CD变化0.000 09; Cm随网格密度的增加呈非单调变化, 从粗网格到密网格, Cm变化约0.000 4.以上计算结果说明, 中等网格已基本消除网格依赖性, 满足本文研究要求.
|
|
下载CSV 表 2 CRM-WBH构型的气动特性(α = 3.0°) Tab.2 Aerodynamic characteristics of CRM-WBH model (α = 3.0°) |
利用CRM-WBHS模型和风洞试验有限元模型, 采用CFD方法和流固耦合计算方法模拟CRM-WBHS模型的绕流流场, 通过与CRM-WBH模型气动特性计算结果和ETW风洞变形测量、测力和测压试验结果的对比分析, 研究支撑装置、静气动弹性变形对CRM翼/身/平尾组合体数值模拟结果的影响.
4.1 计算网格在CRM-WBH模型中等网格的基础上, 构造了CRM-WBHS模型的多块对接结构网格, 局部网格拓扑及表面网格见图 5.为尽量避免网格拓扑和网格分布引起的计算结果差异, 在前机身、翼身结合部、机翼和尾翼附近采用了与CRM-WBH模型相同的网格拓扑及网格分布.整体计算网格单元数达到6.061×107, 物面第1层网格法向无量纲距离y+≈0.67.
|
| 图 5 CRM-WBHS模型网格拓扑和对称面网格(局部) Fig.5 Grid topology and surface grids of CRM-WBHS model (local) |
图 6给出了采用流固耦合计算方法得到的CRM-WBHS模型机翼后缘挠度(dy)和扭转角(dθ)随来流迎角的变化, 同时给出了ETW风洞试验的变形测量结果, 其中η为机翼展向无量纲距离.来流条件为: Ma = 0.85, Re = 5.0×106, α = 0.0~4.0°, 速压q = 60.534 kPa, 载荷因子q/E = 3.342×10-7.在相同迎角下, 机翼后缘的挠度和扭转角由翼根到翼梢逐渐增加; 在相同机翼展向站位, 机翼后缘的挠度和扭转角随来流迎角的增加而逐渐增加.与ETW风洞试验的变形测量结果相比较, α≥3.0°以后, 机翼后缘挠度和扭转角在外翼站位(η = 0.6~1.0)略低于试验结果, 但总体而言, 计算结果与试验结果吻合良好.从以下的分析中可以看出, 机翼静气动弹性变形的上述变化规律将对机翼上的压力分布和气动特性随迎角的变化趋势产生显著影响.
|
| 图 6 机翼后缘的挠度和扭转角沿展向位置分布 Fig.6 Spanwise wing bending and twist distributions at wing trailing edge |
图 7给出了Ma = 0.85, Re = 5.0×106, α = 3.0°时, 采用CFD方法(CRM-WBHS_CFD)和流固耦合计算方法(CRM-WBHS_FSC)得到的CRM-WBHS模型机翼4个典型展向位置的压力系数(Cp)分布曲线, 其中横坐标x采用当地弦长无量纲的流向坐标, 同时还给出了CRM-WBH模型CFD计算结果(CRM-WBH_CFD)和ETW风洞试验的测压结果.
|
| 图 7 CRM翼/身/平尾组合体典型展向站位压力系数分布 Fig.7 Pressure coefficients distributions at different spanwise locations of CRM wing/body/tail model |
对比CRM-WBHS_CFD与CRM-WBH_CFD的计算结果可以看出, 模型支撑装置对压力分布的影响从翼根到翼梢一直存在, 使得机翼上翼面的激波位置前移, 对机翼上翼面其它位置和机翼下翼面的压力分布基本没有影响.从CRM-WBHS_FSC与CRM-WBHS_CFD的对比可以看出, 静气动弹性变形对压力分布的影响由翼根到翼梢逐渐增加.从靠近翼根站位(η = 0.131)直到机翼中部站位(η = 0.502), 静气动弹性变形主要使得机翼上表面激波位置略微提前, 对机翼上表面的其它位置和机翼下表面的压力分布基本没有影响; 在η = 0.846站位, 静气动弹性使得上翼面激波位置以前的负压显著下降、激波位置前移, 下翼面50%弦长以前的压力下降; 在靠近翼梢站位(η = 0.950), 静气动弹性的影响使得机翼上翼面呈现出明显的双激波结构, 并导致上翼面50%弦长以前负压下降、下翼面50%弦长以前的压力下降. CRM-WBHS_FSC压力分布计算结果更加接近ETW风洞测压试验结果.
为进一步考察模型支撑和静气动弹性变形对机翼气动载荷的影响, 图 8给出了α = 3.0°时, CRM翼/身/平尾组合体模型沿机翼展向剖面升力系数(CLs)分布, 同时给出了ETW的试验结果.与CRM-WBH_CFD计算结果相比较, 模型支撑装置导致机翼剖面升力系数略微下降; 模型静气动变形使得剖面升力系数显著下降并更接近试验值.
|
| 图 8 CRM翼/身/平尾组合体展向剖面升力系数分布 Fig.8 Spanwise section lift coefficients distributions of CRM wing/body/tail model |
图 9给出了CRM-WBHS模型CFD计算(CRM-WBHS_CFD)和流固耦合计算(CRM-WBHS_FSC)得到的纵向气动特性曲线, 同时还给出了CRM-WBH模型CFD计算结果(CRM-WBH_CFD)以及ETW风洞试验结果.来流条件为: Ma = 0.85, Re = 5.0×106, α = 0.0~4.0°; 速压q = 60.534 kPa, 载荷因子q/E = 3.342×10-7.
|
| 图 9 CRM翼/身/平尾组合体的气动特性 Fig.9 Aerodynamic characteristics of CRM wing/body/tail model |
对比CRM-WBHS_CFD与CRM-WBH_CFD可以看出, 模型支撑装置使得升力系数、阻力系数下降, 俯仰力矩系数增加, 升力线斜率基本不变, 模型支撑装置对气动特性计算结果的影响量基本不随迎角变化. CRM-WBHS_FSC与CRM-WBHS_CFD的对比可以看出, 在计算迎角范围内, 静气动弹性变形使得升力系数、阻力系数进一步下降, 俯仰力矩系数进一步增加, 升力线斜率略有减少、静稳定性下降; 静气动弹性变形对升力系数和阻力系数的影响量随迎角的增加而增加; α≤3.0°时, 静气动弹性对俯仰力矩系数的影响量随迎角的增加而增大, α>3.0°以后, 影响量随迎角的增加而逐渐减少.采用流固耦合方法得到的CRM-WBHS模型升力系数和阻力系数的计算结果更加接近试验值; 俯仰力矩系数计算结果与试验结果的吻合程度得到进一步改善, 但依然有较大差距.数值计算结果与风洞试验结果之间的差异需从风洞试验数据的处理和数值计算方法两个方面进一步开展研究工作.
5 结论采用CFD方法和流固耦合方法计算了CRM-WBHS构型的气动特性, 通过与CRM-WBH构型计算结果和ETW风洞试验结果的比较, 研究风洞模型支撑装置和静气动弹性变形对CRM翼/身/平尾组合体气动特性的影响, 主要结论如下:
(1) 计算模型中包含支撑装置使得机翼上翼面激波位置前移; 机翼展向气动载荷略微下降; 支撑装置对升力系数、阻力系数和俯仰力矩系数的影响量随来流迎角的增加基本是常量.
(2) 计算模型中包含静气动弹性变形使得内侧机翼上表面激波位置前移, 外侧机翼上表面激波位置前移并导致激波前负压减少, 机翼翼梢上表面激波结构发生变化; 机翼展向气动载荷显著下降; 静气动弹性变形对升力系数和阻力系数的影响量随来流迎角的增加而增加, 对力矩系数的影响量随来流迎角的增加呈非线性变化.
(3) 计算模型中同时包含模型支撑装置和静气动弹性变形, 显著降低了气动特性计算结果与试验结果之间的差异.
致谢 感谢张玉伦、王光学、杨小川等同志在TRIP软件并行计算及网格生成方面所作的研究工作.| [1] |
Vassberg J C, Dehaan M A, Rivers S M, et al. Development of a common research model for applied CFD validation studies[R]. AIAA 2008-6919, 2008.
|
| [2] |
Rivers M B, Dittberner A. Experimental investigations of the NASA common research model in the NASA Langley national transonic facility and NASA Ames 11-Ft transonic wind tunnel (invited)[R]. AIAA 2011-1126, 2011.
|
| [3] |
Rivers M B, Quest J, Rudnik R. Comparison of the NASA common research model European transonic wind tunnel test data to NASA test data (invited)[R]. AIAA 2015-1093, 2015.
|
| [4] |
Vassberg J C, Tinoco E N, Mani M, et al. Summary of the fourth AIAA computational fluid dynamics drag prediction workshop[J]. Journal of Aircraft, 2014, 51(4): 1070-1089. DOI:10.2514/1.C032418 |
| [5] |
Rivers M B, Hunter C A. Support system effects on the NASA common research model[R]. AIAA 2012-0707, 2012.
|
| [6] |
Rivers M B, Hunter C A, Campbell R L. Further investigation of the support system effects and wing twist on the NASA common research model[R]. AIAA 2012-3209, 2012.
|
| [7] |
Keye S, Rudnik R. Validation of wing deformation simulations for the NASA CRM model using fluid-structure interaction computations[R]. AIAA 2015-0619, 2015.
|
| [8] |
王运涛, 王光学, 张玉伦. TRIP2.0软件的确认:DPW Ⅱ复杂组合体的数值模拟[J]. 航空学报, 2008, 29(1): 34-40. Wang Y T, Wang G X, Zhang Y L. Validation of TRIP2.0:numerical simulation of DPW Ⅱ complex con-figuration[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(1): 34-40. DOI:10.3321/j.issn:1000-6893.2008.01.005 (in Chinese) |
| [9] |
王运涛, 李伟, 李松, 等. 梯形翼风洞试验模型数值模拟技术[J]. 航空学报, 2016, 37(4): 1159-1165. Wang Y T, Li W, Li S, et al. Numerical simulation of trapezoidal wing wind tunnel model[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(4): 1159-1165. (in Chinese) |
| [10] |
Van Leer B. Towards the ultimate conservative difference scheme. Ⅱ. Monotonicity and conservation combined in a second-order scheme[J]. Journal of Computational Physics, 1974, 14(4): 361-370. |
| [11] |
Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605. DOI:10.2514/3.12149 |
| [12] |
Yoon S, Jameson A. Lower-upper symmetric-Gauss-Seidel method for the Euler and Navier-stokes equations[J]. AIAA Journal, 1988, 26(9): 1025-1026. DOI:10.2514/3.10007 |
| [13] |
Harder R L, Desmarais R N. Interpolation using surface splines[J]. Journal of Aircraft, 1972, 9(2): 189-191. |
| [14] |
Rendall T C, Allen C B. Efficient mesh motion using radial basis functions with data reduction algorithms[J]. Journal of Computational Physics, 2009, 228(17): 6231-6249. DOI:10.1016/j.jcp.2009.05.013 |
| [15] |
Soni B K. Grid generation for internal flow configura-tions[J]. Computers & Mathematics with Applications, 1992, 24(5/6): 191-201. |
| [16] |
孙岩, 邓小刚, 王运涛, 等. RBF_TFI结构动网格技术在风洞静气动弹性修正中的应用[J]. 工程力学, 2014, 31(10): 228-233. Sun Y, Deng X G, Wang Y T, et al. Application of structural dynamic grid method based on RBF_TFI on wind tunnel static aero-elastic modification[J]. Enginee-ring Mechanics, 2014, 31(10): 228-233. (in Chinese) |