2. 民用飞机设计数字仿真技术北京市重点实验室，北京 102211
2. Beijing Commercial Aircraft Design and Numerical Simulation Laboratory, Beijing 102211, China
自20世纪60年代起, CFD技术在飞机设计过程中扮演了越来越重要的角色, 目前已经成为飞机设计的主要手段, 国外主要飞机制造商都拥有具有完全自主知识产权的核心求解器.例如, 民用航空制造业龙头企业美国波音公司涉及气动专业的工作就是依托TRANAIR软件开展的, 该工具系波音公司自主开发, 应用历史已经长达30多年, 并且仍然在持续开发过程中.相较而言, 我国在这一方面还存在明显差距, 设计工具主要依赖商用软件, 无法实现自主可控, 尽管各科研院所、高校也有自编代码, 但难以全面满足工程需求, 其原因在于工程对工具有较为严苛的要求, 一款CFD工具想要投入工程应用至少须满足以下5方面要求:
(1) 计算快速:计算效率是工程上至关重要的考量指标, 由于飞机设计过程中计算量大、计算类目繁多, 诸如全包线飞行工况、广泛气动设计寻优空间、大量多学科耦合迭代分析等, 工程上往往希望针对全机级CFD分析时长能压缩至分钟级甚至更低, 如颤振分析非定常计算效率最好能压缩至秒级, 当然这对CFD求解器提出了不小挑战.
(2) 分析准确:准确的流场分析对于飞机设计至关重要, 流场分析结果不可信将严重影响设计方案可靠性, 因此CFD工具必须准确捕捉诸如激波、边界层效应等流场特征, 分析结果要通过大量标准模型风洞试验验证.当然, 计算精准度与效率往往存在取舍, 一般而言, 计算网格越密集则分析越准确, 反之亦然.
(3) 操作简单:工程人员希望CFD工具尽可能易于操作. CFD分析前处理工作中至关重要的就是人工划分网格工作, 网格质量优劣直接影响数值求解精度, 而网格划分不仅要求工程人员具备足够经验, 还耗费大量时间、精力, 由于有人工介入也不利于自动化分析.因此工程对从前处理、网格生成、流场求解到后处理全流程自动化的工具有迫切需求.
(4) 运行稳定:鲁棒性是一款工具能够投入工程使用的必要条件, 即要求工具能够适应各类飞行工况, 主要挑战在于失速点附近状态的数值模拟, 由于采用Reynolds时均湍流模型的方法已经不能准确表征强分离物理现象, 所以对于CFD工具而言过失速状态流场计算是一个不小的挑战.
(5) 功能全面:工程要求CFD工具功能齐全, 能够涵盖各类设计状态.从飞机构型上划分包括常规构型、翼身融合构型、高升力构型等; 从流动速度上划分包括亚声速、跨声速、超声速流动等; 从应用领域上划分包括气动布局评估、气动设计、气动特性计算、气弹分析、飞行载荷计算以及多学科耦合设计等.
基于上述需求, 中国商飞北研中心开发了一款面向工程应用的CFD快速自动分析工具SUN程序, 主要针对商用飞机亚声速和跨声速飞行特点, 可对机身、机翼、挂架和短舱部件组合的任意构型飞机流场进行快速分析, 代码行数约4×104行, 采用准面向对象程序设计架构和OpenMP多线程并行技术, 实现结构化贴体网格生成、网格投影、网格嵌套、无黏外流计算、边界层耦合和后处理全流程一体化分析, 具有自动化程度高、建模简单、计算快速、结果准确和鲁棒性强等特点, 对于百万量级的定常流场计算, 在单机上计算时长约6 min, 适用于飞机气动布局评估、气动设计、气动特性计算、气弹分析、飞行载荷计算以及多学科耦合设计等相关工作, 目前已经在中国商飞公司多个项目中得到应用.
1 程序架构SUN程序采用Fortran语言进行开发并利用module模块实现类似于C++语言的面向对象程序设计, 图 1展示了SUN程序具体架构图, 该模式有利于程序开发团队共同作业与维护. SUN程序主要模块包括数学计算、几何前处理、网格生成、网格嵌套、流场求解、流动后处理6大模块, 每个模块包含1个主调用程序(以Main_+模块名字段开头)以及Lib和Src两个文件夹目录, Lib文件夹存放程序所需要调用的类文件, 实现对象接口统一管理, Src目录存放库文件具体调用的源码文件, 每段代码文件行数保持在200行左右, 各子模块通过Makefile文件编译打包成库, 主程序链接这些库完成最终编译, 程序提供Windows/Linux/Mac OS操作系统3种版本.
|
| 图 1 SUN程序架构图 Fig.1 Framework of SUN code |
相比于非结构化网格, 结构化贴体网格[3]具有网格节点相邻关系明确、隐式或多阶格式内存消耗小、附面层流动易于模拟等优势.但针对飞机构型复杂的几何形态, 很难用统一的方式自动生成高质量结构化贴体网格, 因此需要采用嵌套网格技术[4-5], 即将整个流场分区, 每个子区根据部件的几何特征独立生成高质量计算网格, 子区网格的边界面建立不受其他子区网格影响, 各区块相对独立进行流场分析并通过重叠区域流场数据传递计入相互间影响, 但在两种几何部件连接位置, 建立的网格存在间隙, 从属部件流场无法捕获主部件物面边界影响, 因此需要把从属部件网格投影到主部件几何上去建立Collar网格[6], 网格自动投影技术系本文作者原创性提出, 在国际范围内尚属首次.
2.1 飞机构型几何描述SUN程序各部件几何物面采用一簇离散点描述, 机身、机翼、挂架、短舱几何物面分别由沿机头到机尾连续截面框离散点、沿展向翼型截面离散点、沿法向翼型截面离散点、沿周向翼型截面离散点表征(见图 2).
|
| 图 2 飞机构型几何描述 Fig.2 Geometric representation of aircraft configuration |
图 3为网格生成模块技术流程图, SUN程序采用Taylor样条曲线拟合[7]的参数化方法描述飞机几何外形, 机身网格采用椭球和圆柱正交网格生成坐标系, 通过外推表面网格得到三维空间计算网格, 机翼/短舱/挂架网格采用C-H网格拓扑结构保角变换生成二维翼型网格, 沿展向线性插值组合出机翼三维网格, 沿法向线性插值组合出挂架三维网格, 沿周向旋转一周生成短舱三维网格.
|
| 图 3 SUN程序网格生成模块概况图 Fig.3 Flow chart of grid auto-generating module |
图 4为SUN程序网格嵌套模块技术[8]流程图, 由于各部件网格独立生成, 因此不同块网格间势必发生交叉.进入到物面内部的点称作几何内点; 与几何内点相邻的点为网格交叉点; 六面体网格单元任一表面只要存在一个常规网格点即为常规表面; 递进的只要存在一个网格交叉点则为交叉表面; 全部为几何内点则该面为几何内表面; 六面体单元只要存在常规网格点即为常规单元; 递进的只要存在网格交叉点即为交叉单元; 全部是几何内点则为几何内单元.由于几何内单元没有物理意义因此需要剔除, 即“挖洞”处理.流场信息以交叉表面对流通量的形式进行交互传递.
|
| 图 4 SUN程序嵌套网格模块概况图 Fig.4 Flow chart of chimera grid auto-checking module |
为了准确捕获主部件物面对从属部件流场的影响, 如机身物面对机翼流场、机翼/短舱物面对挂架流场、机身/短舱物面对尾吊布局飞机支撑件影响等, SUN程序创新提出网格自动投影方法, 将部件几何连接区域从属部件网格(如机翼、垂尾、平尾、立尾等)投影到主部件表面(如机身、翼身融合布局主升力体)以划分Collar网格, 图 5和图 6分别展示了SUN程序机翼网格自动投影到机身物面[9]、挂架网格自动投影到机翼/短舱物面[10]的流程图, 这两项技术均已申报国家发明专利.
|
| 图 5 SUN程序机翼网格投影到机身物面流程图 Fig.5 Flow chart of wing mesh auto-projecting module |
|
| 图 6 SUN程序挂架网格投影到机翼/短舱物面流程图 Fig.6 Flow chart of pylon mesh auto-projecting module |
为了满足工程对计算效率要求, SUN程序采用耦合外侧无黏流动和黏性边界层流动方法[11-13]开展流场模块开发, 具体技术流程见图 7, 该方法能够准确模拟激波和边界层效应.由于外流网格无需与边界层网格保持同等密度, 相比RANS求解网格总量级低1~2阶, 计算效率高且只需单机运行, 无需高性能计算集群的硬件支持, 计算效率优势明显, 全机构型在粗网格和密网格上流场分析时长分别可控制在1 min和5 min左右.
|
| 图 7 SUN程序流场求解模块概况图 Fig.7 Flow chart of flow simulating module |
对于跨声速流场分析, 需要准确捕获激波, 因此外侧无黏绕流采用守恒形式Euler方程描述, 该方法相较面元法、速势方程法能够准确处理激波出现后的流动非线性问题, 数值离散采用格心格式有限体积法, 引入隐式迎风格式[14]提高收敛判据, 采用二阶Roe格式[15]计算Riemann问题对流通量, 采用LU+GMRES方法[16]和OpenMP并行技术加速流场更新, 详细理论推导可参见文献[17].
3.2 内侧边界层流动边界层方程[18]数值离散采用格点格式有限差分法, 应用Cebeci-Smith零方程模型[19]模拟湍流平均流动, 对于强分离流动区域采用半反模式内外流耦合技术[20-21], 通过Keller Box方法[22]离散边界层前缘驻线方程和通用方程并采用Newton方法求解, 对于尾迹区域流动采用Green积分边界层方法[23], 详细理论推导参见文献[24].
3.3 内外流耦合技术由于飞机在真实飞行过程中Reynolds数非常高, 外侧绕流中的黏性效应可以忽略不计, 但边界层内部黏性剪切效应将导致流场质量流溢出, 因此须把溢出效应修正到外侧无黏绕流, 文中采用可穿透边界条件模拟边界层溢出效应.
3.4 静气动弹性分析模块以SUN程序为核心求解器搭建基于CFD/CSD松耦合的静气动弹性分析平台, 先采用径向基函数RBF方法将气动载荷插值到结构有限元模型, 再采用商用NASTRAN软件进行静力变形分析, 然后再将结构变形传递回气动模型, 循环迭代直至静气弹变形收敛到容许值, 图 8展示了静气动弹性分析流程.
|
| 图 8 SUN程序静气动弹性分析解模块流程图 Fig.8 Flow chart of static aero-elastic analyzing module |
本节将展示SUN程序对一些标准模型的分析结果以及与风洞试验数据的对比验证, 包括NACA0012翼型，DLR F6标准模型飞机，NASA CRM标准模型飞机.
4.1 NACA0012标准模型图 9展示了不同工况下采用SUN程序对标准模型NACA0012翼型的计算和风洞试验对比验证结果, 飞行工况包括: (1)Mach数0.3，攻角3.54°，Reynolds数1.68×106; (2)Mach数0.5，攻角5.23°，Reynolds数2.91×106; (3)Mach数0.749，攻角2.00°, Reynolds数1.17×107.可以看到, 3组数据与风洞试验结果都有较好的贴合度, 最后一组数据在激波后位置处与风洞试验稍有差别, 这是由于文中采用的Cebeci-Smith零方程湍流模型不能较好模拟激波诱导后的流动分离, 分离泡跟试验相比偏小, 而分离泡会导致流动加速, 压力值减小, 因此数值模拟的加速效应较试验偏弱.
|
| 图 9 NACA0012标准模型流场分析对比验证 Fig.9 NACA0012 model flow simulation verification |
图 10展示了Mach数0.5，Reynolds数2.9×106工况下采用SUN程序计算的NACA0012翼型升阻力曲线与风洞试验结果对比验证.可以看到, 升力曲线基本与试验数据吻合, 阻力曲线由于数值小、误差大, 与试验值稍有偏差, 但基本上满足工程应用要求.
|
| 图 10 NACA0012标准模型升阻力曲线验证 Fig.10 NACA0012 model lift and drag curve verification |
图 11展示了飞行工况Mach数0.75, 飞行攻角1.23°，Reynolds数3×105，应用SUN程序对DLR F6标准模型飞机流场分析以及与风洞试验数据对比验证结果, 流场计算基本与风洞试验数据吻合, 本算例采用双线程并行在Intel(R) Xeon(R) CPU E5-2667单机上运行耗时约4 min 32 s.
|
| 图 11 DLR F6标准模型流场分析对比验证 Fig.11 DLR F6 model flow simulation verification |
图 12展示了飞行工况Mach数0.85, 飞行攻角2.013 29°, Reynolds数3×107, 参考面积191.553 5 m2应用SUN程序对NASA CRM标准模型飞机[25]流场分析、静气动弹性分析以及与风洞试验数据对比验证结果, 本算例采用3线程并行单机运行耗时约5 min 40 s, 静气动弹性耦合迭代分析耗时约40 min.图中红色三角表示风洞试验数据, 黑色曲线表示按照原始模型进行流场分析后计算数据, 可以看到越靠近外翼段其与风洞试验数据偏差越大, 这是由静气动弹性效应引起的.蓝色曲线表示对原始模型进行变形修正后的CFD分析结果, 紫色曲线表示直接对原始模型进行静气动弹性分析后结果, 这两组数据与风洞试验结果基本吻合, 验证了SUN程序流场计算和静气动弹性分析的准确度.
|
| 图 12 NASA CRM标准模型流场分析对比验证 Fig.12 NASA CRM flow simulation verification |
(1) SUN程序流场计算快速、准确、可靠、鲁棒性强、工程应用前景广阔.
(2) 基于无黏/黏性边界层耦合流场计算方法能够准确模拟激波和边界层效应, 单机计算效率优势明显, 但尚不能分析较强分离流动, 系下一步研究重点.
(3) 嵌套网格区域流场守恒性是亟待解决的难点, 目前的不守恒算法导致不同网格参数下力系数发生微幅变动.
(4) 对于采用嵌套网格的流场求解方法必须无限加密几何连接区域网格(计算效率降低)或建立Collar网格, 否则流场分析结果精准度差.
(5) 静气动弹性分析误差来源包括CFD分析误差、CFD/CSD载荷插值误差、静力分析误差、CSD/CFD位移插值误差、变形误差等, 与风洞试验完全吻合尚存在一定难度.
| [1] |
Samant S S, Bussoletti J E, Johnson F T, et al. TRANAIR-a computer code for transonic analyses of arbitrary configurations[R]. AIAA 1987-0034, 1987.
|
| [2] |
孙明哲, 李政德, 李丽雅, 等.飞机定常流场快速自动分析SUN程序软件著作权(申请中). Sun M Z, Li Z D, Li L Y, et al. Fast and automatic steady flow solver SUN software copyright(in Chinese). |
| [3] |
Thompson J F, Warsi Z U, Mastin C W. Numerical Grid Generation[M]. New York: North Holland, 1985.
|
| [4] |
Steger J L, Dougherty F C, Benek J A. A chimera grid scheme[C]. Ghia K N, Ghia U. Advances in Grid Generation. New York: ASME, 1983: 59.
|
| [5] |
Benek J A, Buning P G, Steger J L. A 3D chimera grid embedding technique[R]. AIAA 1985-1523, 1985.
|
| [6] |
Parks S J, Buning P G, Steger J L, et al. Collar grids for intersecting geometric components within the chimera overlapped grid scheme[R]. AIAA 1991-1587, 1991.
|
| [7] |
郭云, 吴松强, 李建蜀. 利用泰勒公式实现样条曲线的拟合[J]. 贵州工业大学学报, 1997, 26(5): 19-22. Guo Y, Wu S Q, Li J S. The algorithm for curve fitting with cubic spline and implement[J]. Journal of Guizhou University of Technology, 1997, 26(5): 19-22. (in Chinese) |
| [8] |
孙明哲.一种嵌套网格建立方法及装置: 201710207741.0[P]. 2017-07-21.有权. Sun M Z. A Chimera Grid Generation Method and Tool: 201710207741.0[P]. 2017-07-21(in Chinese). |
| [9] |
孙明哲, 李政德.一种网格投影方法及装置: 201710648144.1[P]. 2017-08-01. Sun M Z, Li Z D. An automatic mesh projection method and tool: 201710648144.1[P]. 2017-08-01(in Chinese). |
| [10] |
孙明哲, 李政德.飞机挂架部件结构化贴体网格自动生成及投影技术[P].专利申请中. Sun M Z, Li Z D. Automated generation and projection of structured body-fitted grids for aircraft pylon components[P]. in application(in Chinese). |
| [11] |
Lighthill M J. On displacement thickness[J]. Journal of Fluid Mechanics, 1958, 4: 383-392. DOI:10.1017/S0022112058000525 |
| [12] |
Cebeci T, Kaups K, Ramsey J A. A general method for calculating three-dimensional compressible laminar and turbulent boundary layers on arbitrary wings[R]. NASA CR-2777, 1977.
|
| [13] |
Williams B R. The prediction of separated flow using a viscous-inviscid interaction method[J]. Aeronautical Journal, 1985, 89(885): 185-197. |
| [14] |
Jameson A, Yoon S. LU implicit schemes with multiple grids for the Euler equations[R]. AIAA 86-0105, 1986.
|
| [15] |
Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computa-tional Physics, 1981, 43(2): 357-372. |
| [16] |
Saad Y, Schultz M H. GMRES:a generalized minimal residual algorithm for solving nonsymmetric linear system[J]. SIAM Journal on Scientific and Statistical Compu-ting, 1986, 7(3): 856-869. DOI:10.1137/0907058 |
| [17] |
孙明哲, 李政德, 李丽雅.翼身飞机嵌套网格生成及跨音速流场快速分析[C].首届全国空气动力学大会, 2018. Sun M Z, Li Z D, Li L Y. Chimera grid auto-generating and transonic flow simulation for wing-body aircraft configuration[C]. First Chinese Aerodynamics Conferen-ce, 2018(in Chinese). |
| [18] |
Howarth L. The boundary layer in three dimensional flow.-part Ⅰ. Derivation of the equations for flow along a general Curved surface[J]. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 1951, 42(326): 239-243. DOI:10.1080/14786445108561259 |
| [19] |
Cebeci T, Smith A M O. Analysis of turbulent boundary layers[M]. New York: Academic Press, 1974.
|
| [20] |
Cebeci T. An inverse boundary-layer method for compressible laminar and turbulent boundary layers[J]. Journal of Aircraft, 1976, 13(9): 709-717. DOI:10.2514/3.44552 |
| [21] |
Le Balleur J C. Strong matching method for computing transonic viscous flows including wakes and separations[Z]. La Recherche Aerospatiale, 1981.
|
| [22] |
Keller H B. Numerical methods in boundary layer theory[J]. Annual Review of Fluid Mechanics, 2003, 10(1): 417-433. |
| [23] |
Green J E. Application of head's entrainment method to the prediction of turbulent boundary layers and wakes in compressible flow[R]. RAE TR 72079, 1976.
|
| [24] |
孙明哲, 林榕婷, 李政德.一种无粘流动与边界层耦合流场分析方法[C]. 2019全国力学大会, 2019. Sun M Z, Lin R T, Li Z D. Fast flow analyzing method based on inviscid-viscous coupling method[C]. Chinese Congress of Theoretical and Applied Mechanics, 2019(in Chinese). |
| [25] |
Vassberg C J, Dehaan M A, Rivers S M, et al. Development of a common research model for applied CFD validation studies[R]. AIAA 2008-6919, 2008.
|