超大规模气动弹性数值模拟软件研制(2017)

王运涛, 孟德虹, 孙 岩, 洪俊武

(1. 中国空气动力研究与发展中心 计算空气动力研究所, 四川 绵阳 621000;2. 中国空气动力研究与发展中心 空气动力学国家重点实验室, 四川 绵阳 621000)

0 引 言

高性能计算的峰值性能已于2008年跨越P级(PetaFlops，每秒千万亿次浮点运算)门槛，目前正在向E级(ExaFlops，每秒百亿亿次浮点运算)进发。自2010年11月以来，以“太湖之光”(峰值性能125.4P，持续性能93.0P)、“天河二号”(峰值性能54.9P，持续性能33.9P)为代表的超级计算机, 已经多次在全球TOP 500超级计算机排行榜名列第一。截止2017年11月，中国超级计算机TOP 500的上榜总数达到历史最高的202台。随着国内超级计算机硬件技术的飞速发展，国内超大规模的计算应用也取得了突破性进展。2016年，中国科学院软件研究所杨超等人[1]在美国盐湖城举行的全球超级计算大会上获得了国际高性能计算应用领域最高奖“戈登贝尔”奖，获奖应用名称为“千万核可扩展全球大气动力学全隐式模拟”。2017年，清华大学付昊桓等人在美国丹佛举行的全球超级计算大会上再次获得“戈登贝尔”奖，获奖应用项目名称为“非线性地震模拟”。随着我国航空航天飞行器的自主创新发展，国内飞行器设计与研制单位正在逐步装备百万亿次量级的高性能计算机，其中中国空气动力研究与发展中心装备的高性能计算机已经达到1.5P。以高保真气动数据库的高效生产和关键气动技术攻关为牵引的计算空气动力学(CFD)正在成为国产高性能计算机持续发展的原始驱动力之一。

高性能计算机的发展是飞行器多学科耦合设计的硬件基础，而多学科耦合的数值模拟工具又是多学科耦合设计的软件基础，气动弹性问题的数值模拟属于典型的多学科耦合数值模拟问题。始于2010年、由AIAA组织的DPW(Drag Prediction Workshop)系列会议已经成功举办了六届，其主要目的是评估现代CFD技术模拟运输机高速构型的阻力预测能力，为CFD技术下一步的发展提出建议和意见。为考核现有流固耦合软件的数值模拟能力，2016年6月召开的第六届DPW组委会将CRM(Common Research Model)翼身组合体模型的流固耦合数值模拟(Case 5)作为可选工况。来自世界各地的25个研究机构采用25种软件提供了48组网格收敛性(Case 1)的数值模拟结果[2]，其中只有4家研究机构提供了4组流固耦合的计算结果。由此可见相对于常规的CFD数值模拟，静气动弹性数值模拟技术尚未成熟[3]。始于2012年的由AIAA组织的AePW(Aeroelastic Prediction Workshop)系列会议已经举办了两届，其主要目的是评估气动弹性数值模拟技术现状，为气动弹性数值模拟技术和相应风洞试验技术的发展提出意见和建议。2016年1月召开的第二届AePW组委会将BSCW(Benchmark Supercritical Wing)模型的颤振模拟(Case2)作为必选工况[4]。来自世界各地的7个研究团队提供了13组流固耦合的数值模拟结果，数值模拟得到的颤振起始动压与试验结果最大相差22%。这说明，颤振数值模拟技术与风洞试验技术本身，均需要进一步开展研究工作。

由中国空气动力研究与发展中心(CARDC)与大连理工大学、北京航空航天大学、中国商飞上海飞机设计研究院、成都飞机设计研究所等五家单位联合申请的国家重点研发计划项目“数值飞行器原型系统开发”于2016年8月正式立项。其中，研制具有完全自主知识产权的超大规模气动弹性数值模拟软件是项目的重要研究内容之一。该软件的研究目标是：基于本项目开发的飞行器空气动力学模拟软件和结构动力学分析软件，突破百亿量级超大规模流场计算网格高效鲁棒变形、海量耦合数据高效高精度传递等技术瓶颈，实现百亿量级网格规模、60万处理器核的超大规模静/动气动弹性数值模拟，构建飞行器气动弹性机理研究和工程应用的软件平台，提高飞行器气动弹性问题数值模拟的效率和预测精度，满足型号研制对气动弹性模拟能力的需求。

本文介绍了超大规模气动弹性数值模拟软件的研制情况，主要包括超大规模结构网格生成技术、基于结构网格的动网格技术、气动弹性数值模拟软件主要功能模块的并行效率和并行一致性测试、软件的确认工作等内容。

1 超大规模气动弹性数值模拟软件简介

超大规模气动弹性数值模拟软件主要包括CFD求解模块、CSM(Computational Structural Mechanics)求解模块、耦合界面数据传递和动态网格变形4个主要功能模块，通过主控程序有序组织上述功能模块，从而实现复杂飞行器气动弹性问题的数值模拟。其中，超大规模并行计算技术贯穿于软件的各个功能模块之中。目前，“数值飞行器原型系统开发”项目组已经开发了超大规模气动弹性数值模拟软件的初级版本，各个功能模块的主要算法简单介绍如下。

(1) CFD求解模块。采用CARDC自主研发的亚跨超CFD软件平台(Trisonic Platform，TRIP)——TRIP软件,经过了系统的验证和确认工作[5-6]并广泛应用于多种飞行器的气动设计与评估。该软件采用有限体积方法离散任意坐标下的雷诺平均Navier-Stokes方程(RANS)，RANS方程无黏项离散采用二阶精度ROE迎风格式[7]，黏性项离散采用二阶中心格式，湍流模型包含Spalart-Allmaras SA一方程模型[8]、Menter’s SST两方程模型[9]，定常问题离散方程组的求解采用LU-SGS方法[10]，非定常问题离散方程组的求解采用Jameson“双时间步”方法[11]，采用低速预处理技术、多重网格技术和基于MPI或MPI/OpenMP的大规模并行技术加速收敛。

(2) CSM计算模块。静气动弹性数值模拟采用柔度矩阵方法求解线性结构静力学方程，获得气动载荷作用下的物面变形；颤振等动气动弹性数值模拟采用模态法求解结构动力学方程。目前，结构柔度矩阵以及模态的计算均基于飞行器有限元模型采用商业结构动力学分析软件获得，适用于飞行器结构线性变化范围内气动弹性问题的模拟，对于结构非线性气动弹性问题的模拟需要进一步开展研究工作。

(3) 耦合界面数据传递模块。采用无限板样条IPS(Infinite Plate Spline)方法或薄板样条TPS(Thin Plate Spline)插值方法[12]构建CFD计算模块与CSM计算模块之间的气动载荷与结构变形传递矩阵。

(4) 动态网格变形模块。CFD计算网格的变形采用基于径向基函数RBF(Radial Basis Functions)[13]与超限插值TFI(Transfinite Interpolation)[14]的复合型动态网格变形方法RBF_TFI[15]。

(5) 静气动弹性计算流程。采用松耦合方式建立静气动弹性模拟的一般流程(图1)。第一步，CFD模块计算出第n迭代步流场Un后，将壁面压力Pn传递给CSM模块；第二步，CSM模块利用输入的气动载荷和边界约束条件，计算出第n+1迭代步的结构变形Wn+1；第三步，将结构变形Wn+1传递给CFD模块；第四步，CFD模块利用变形后外形，更新计算网格，计算第n+1迭代步的流场Un+1；然后循环上述过程，不断得到下一迭代步的位移和流场，直至变形位移和流场均达到收敛。

图1 静气动弹性问题模拟流程Fig.1 Outline of static aero-elastic problem simulation

(6) 动气动弹性计算流程。采用松耦合方式建立动气动弹性模拟的一般流程(图2)。第一步，CFD模块计算出t0时间步流场Ut0后，将壁面压力Pt0传递给CSM模块；第二步，CSM模块利用输入的气动载荷和边界约束条件，基于模态迭加的结构动力学方程，采用四步Runge-Kutta方法计算出t0+Δt时间步的结构变形Wt0+Δt；第三步，将结构变形Wt0+Δt传递给CFD模块；第四步，CFD模块利用变形后外形，更新流体计算网格，采用“双时间步”方法计算t0+Δt时间步的流场Ut0+Δt。循环上述过程，不断得到下一时刻的位移和流场，直至计算过程满足收敛判据。由于CFD求解模块和CSM求解模块在数据交换时存在一个时间步的延迟，松耦合方式在时间方向只具有一阶时间精度。下一步将采用预估-校正技术或将CSM求解模块同样构造成含有子迭代的求解方式以提高动气动弹性问题数值模拟的时间精度。

图2 动气动弹性计算流程Fig.2 Outline of dynamic aero-elastic problem simulation

2 超大规模气动弹性数值模拟软件研制

2.1 超大规模静态结构网格生成技术

超大规模静态网格生成技术是实现复杂构型超大规模气动弹性问题数值模拟的前提。以包含机翼/机身/挂架/短舱/平尾/立尾的典型运输机巡航构型为例，目前模拟上述复杂构型的结构网格规模一般在千万量级，在配置较高的微机或工作站上利用商业网格生成软件即可完成。而构造数亿量级乃至百亿量级的结构网格已经远远超出了当前单机能力范围，必须采用新的网格生成策略。超大规模静态网格生成的关键技术主要包括物面网格的高保真度技术，高质量、鲁棒的空间网格构造技术等方面。

超大规模结构网格的生成步骤是：1) 采用商业软件构造千万量级高质量初始结构网格；2) 利用初始结构网格拓扑分块信息，建立网格窗口的信息数据库；3)基于网格窗口的信息数据库创建不同网格窗口之间的关联关系；4) 根据网格规模需求，对棱线网格点进行加密，重新计算棱线网格点分布及棱线网格点坐标，通过超限插值(TFI)生成加密后的窗口网格；5) 采用线性投影方法，对边界属性标记为物面的窗口网格重新向物面投影，获得加密后的物面网格；6) 利用三维TFI生成加密后的体网格。上述超大规模网格生成策略与待模拟外形的复杂程度无关，所构造的超大规模结构网格拓扑结构和网格单元质量与初始网格完全相同，并可以对初始网格进行不同程度的粗化或细化。采用上述方法，可以快速获得百亿量级以上的高质量结构网格，同时可以构造一系列的不同规模网格开展网格收敛性研究。

按照前述方法，对CRM-WBH构型进行了验证性网格重构。初始网格规模为9 237万网格单元，采用网格粗化和细化的方法获得了105 219万(10亿量级)、591 187万(60亿量级)和1 154 662万(115亿量级)网格单元的不同规模网格。图1给出了115亿网格单元的CRM-WBH构型空间网格拓扑。表1给出了不同规模网格的主要信息，其中Ratio表示棱线网格单元加密倍数，Nnode表示网格结点总数，NBL表示边界层网格的法向网格点数目，y+表征第一层网格法向无量纲距离，Nblock表示总的网格块数目，Size表示网格文件所占存储量。

图3 CRM-WBH构型的网格拓扑(115亿)Fig.3 Grid topology of CRM-WBH configuration (11.5 billion)

ParameterRatio=1.0Ratio=2.25Ratio=4.0Ratio=5.0Nnode92 372 9921 052 186 1125 911 871 48811 546 624 000NBL97217385481y+0.50.220.1250.1Nblock73737373Size/MB2 27725 556142 980278 829

2.2 超大规模结构网格变形技术

高质量的网格变形技术是CFD领域的研究热点，也是气动弹性数值模拟的关键技术之一。高质量、鲁棒高效的网格变形技术应满足以下要求：1)自动化。减少甚至避免网格变形过程中的人工干预是气动弹性问题高效数值模拟的基本要求。2)并行化。网格变形过程的并行化是提高网格变形效率有效途径，对于超大规模网格的动气动弹性问题数值模拟，网格变形并行化的需求尤其迫切。3) 适用性。网格变形方法要能够处理各种复杂外形、大变形，以及各种复杂变形(包括不规则变形)。4)通用性。网格变形方法不仅适用于结构网格，还适用于各种类型的非结构网格及混合网格。

本文选择了基于径向基函数的超限插值方法(RBF_TFI)实现复杂构型结构网格的自动变形，利用删减基点后的径向基函数插值更新结构网格块的棱线坐标，然后利用TFI插值重新生成计算网格。该方法的具体实施步骤是：1)根据输入的多块结构网格信息，在物面网格处选择径向基函数基点，构造可逆的插值矩阵；2)根据输入的物面变形确定物面径向基基点位移，利用可逆的插值矩阵确定径向基插值系数；3)采用径向基插值方法更新物面与空间网格的网格棱线坐标；4)采用二维、三维TFI方法构造各个网格块的面网格与体网格。RBF_TFI方法凭借优异的性能在很多程序和算例中得到了广泛应用，但在一些比较极端的情况下，边界层网格第一层高度十分小的时候(这在高雷诺数流动问题和超大规模网格问题中经常遇到)，RBF_TFI容易引起边界层区域网格棱线的扭曲现象，并造成物面附近网格的交叉破坏问题，从而导致流场计算无法继续进行。基于不同算例对该现象进行了深入的分析，认为是物面网格点与空间网格点采用不同的坐标更新方法导致的变形不一致引起了网格交错破坏。通过现象、原因等方面的详细分析，提出一种解决该问题的改进方法，改进后的网格更新流程见图4。

图4 RBF_TFI结构网格变形工作流程Fig.4 Workflow of structured grid deformation with RBF-TFI method

在颤振问题的数值模拟中，为了保证动气动弹性问题数值模拟时间上的二阶精度，流体力学方程求解和结构力学方程求解一般采用真实时间步加上子迭代的方式，在每一个子迭代步均需更新流场计算的网格。在这种情形下，CFD计算网格的更新效率至关重要，尤其是对于亿量级或数十亿量级网格的流场计算。本文采用基于MPI的并行处理技术，开展了结构网格变形并行化技术研究，初步建立了结构网格并行变形的一般流程。具体实现方法如下：首先主进程收集所有进程物面变形，计算全局插值系数，然后将插值系数播撒到所有子进程，每个核利用插值系数矩阵独立采用RBF或RBF+TFI方法更新体网格。物面网格的更新方法与体网格的更新方法相同。

2.3 超大规模气动/结构数据并行传递技术

气动弹性问题数值模拟时，需要将CSM模块获得气动载荷作用下的结构变形传递给动态网格处理模块，进而开展下一个迭代步的CFD计算；将CFD模块获得的气动载荷传递给CSM模块以开展结构分析，位移数据的传递与气动载荷数据的传递均在物面进行。这涉及两个方面的问题，一个是数据传递的精度与守恒性问题，另一个是针对超大规模气动特性问题模拟的效率问题。第一个问题已经在前期工作中得以解决[15-16]。针对第二个问题，为了使得数据传递的效率与气动特性问题数值模拟的效率相匹配，采用了基于MPI的并行处理技术，并行化的耦合数据传递基于径向基函数插值技术实现，与动态变形网格并行化的思路十分相似。具体实现方法是：首先在主结点实现插值系数方程组的求解，获得全局耦合插值的系数，然后将插值系数播撒到各个CPU核，每个核独立实现耦合变量的传递，最后组合实现全局耦合界面数据的高效传递。

3 气动弹性数值模拟软件的并行测试

2017年度，采用AGARD445.6机翼对超大规模气动弹性数值模拟的动网格模块开展了并行效率测试；采用CRM-WBH模型，对CFD求解模块开展了并行一致性测试及并行效率测试；采用AGARD445.6机翼，对动气动弹性模拟模块开展了并行效率测试。

3.1 CRM-WBH模型与AGARD445.6颤振模型

CRM-WBH模型是2010年召开的第四届AIAA阻力预测研讨会(DPW IV)的基准研究模型[17]，设计马赫数为0.85，设计升力系数为0.50。该模型先后在NASA Langley的NTF风洞和NASA Ames的TWT风洞中开展了多种构型的风洞试验。风洞试验模型采用安装于机身后体的叶片尾撑方式固定于风洞迎角变换装置，试验结果包括了气动特性、表面压力分布等。DPW IV会议上数值模拟得到的气动特性与相应风洞试验结果的对比并不理想[18-19]，主要表现：相同迎角下，数值模拟得到的升力系数普遍大于试验结果；相同升力系数下，数值模拟得到的低头力矩系数普遍大于试验结果。

AGARD445.6机翼是用于考核跨声速颤振计算方法和软件的标准模型[20]，风洞试验是在美国NASA Langley 研究中心的跨声速动态风洞完成的。AGARD445.6机翼弦向翼剖面为NASA65A004翼型，试验模型以材质均匀的桃花心木薄板制成。机翼根弦长0.5587 m、展长为0.762 m，机翼展弦比1.6525，根梢比1.5207，1/4弦长后掠角为45°。

3.2 动网格模块并行效率测试

表2给出了AGARD445.6机翼2.05亿结构网格动网格模块并行效率测试结果。采用128～4096个处理器测试了动网格模块的并行效率，网格更新时间相对于流场单步子迭代的计算时间控制在14%以内。在广州超算中心的天河二号集群上开展了十亿量级网格变形并行效率测试。图5给出了AGARD445.6机翼模型16亿多块结构网格变形并行测试的结果，图中横坐标为采用的处理器核数，纵坐标为网格更新时间与单步CFD流场计算时间的比值。对于16亿的流场计算网格，从2 000核到25 000核，网格变形更新时间占单次CFD计算时间的比例基本稳定在15%左右。上述测试表明，采用的动网格并行策略是可行的，并行效率不仅满足工程问题的效率需求，也为数十亿网格规模的气动弹性问题数值模拟提供了一个可能的解决方案。

表2 2.05亿结构网格变形并行效率测试Table 2 Parallel efficiency test of 20.5 million structured grid deformation

图5 16亿结构网格变形并行效率测试Fig.5 Parallel efficiency test of 1.6 billion structured grid deformation

3.3 静气动弹性模块并行效率及一致性测试

采用CRM-WBH构型、1.8亿网格，开展了CFD模块的并行效率测试。表3给出了采用单重网格、960个处理器核到9600个处理器核并行计算相对于96个处理器核的加速比和并行计算效率。在同构MPI并行模式下，9600个处理器核的并行效率达到70%左右。并行一致性是大规模并行计算的另一个关键问题，采用CRM-WBH构型开展了CFD软件的并行一致性测试。计算来流条件为：Ma=0.85，Re=5.0×106，α=3.0°。表4给出了96～9600处理器核并行计算得到的CRM-WBH构型的气动特性。当处理器规模由96个增加到9600个时，升力系数的变化量为0.12%、阻力系数的变化量为0.34%、俯仰力矩系数的变化量为1.1%。计算结果一致性良好。

表3 1.8亿结构网格CFD并行效率测试Table 3 Parallel efficiency test of CFD with 18.0 million structured grid

表4 1.8亿结构网格CRM-WBH构型并行一致性测试Table 4 Parallel consistency test of CRM-WBH configuration with 18.0 million structured grid

3.4 动气动弹性模块并行效率测试

采用AGARD445.6机翼16亿规模的多块对接结构网格，在广州超算中心天河二号超级集群上进行了动气动弹性模块的大规模并行效率测试，并行规模为1536～24 576个处理器核。图6给出了并行效率和加速比，由图可以看出，采用MPI并行方式，相对于1536个处理器核，并行规模为24 576个处理器核时，并行效率达到了77.63%。

图6 AGARD445.6机翼颤振计算并行效率Fig.6 Parallel efficiency of flutter simulation for AGARD 445.6 wing

4 气动弹性数值模拟软件的确认

在前期工作中，开展了HIRENASD机翼模型、DLR-F6模型[16]，CRM-WB模型[21]静气动弹性模拟测试，初步确认了气动弹性模拟软件的静气动弹性模拟精度。本文采用CRM-WBH构型进一步开展了静气动弹性数值模拟精度确认，采用AGARD445.6机翼开展了动气动弹性数值模拟精度确认。

4.1 CRM-WBH构型静气动弹性数值模拟

采用CRM-WBH模型、支撑装置数模及风洞试验模型的有限元模型开展了CRM-WBH构型的静气动弹性模拟。以下为叙述方便，将包含尾撑装置的CRM翼/身/平尾组合体构型简称为CRM-WBHS。图7给出了包含支撑装置的CRM-WBHS构型多块对接结构网格拓扑及表面网格，网格规模6061万。图8给出了CRM-WBH风洞试验模型的有限元模型。

图7 CRM-WBHS构型网格拓扑及表面网格Fig.7 Grid topology and surface grid for CRM-WBHS configuration

图8 CRM-WBH风洞试验有限元模型Fig.8 Finite element model CRM-WBH wind tunnel test

图9为不同迎角下静弹性计算最大结构变形随迭代步数的变化曲线，不同迎角下的计算结果都可以看出,在迭代到第5步后，最大变形量基本不变。

图10给出了CRM-WBH构型气动特性随迎角的变化，计算来流状态：Ma=0.85，Re=5.0×106，α=0°～4.00°；流固耦合计算时，速压q=61.062kPa，载荷因子q/E=3.342×10-7,E为弹性模量。其中，CRM-WBHS_CFD代表采用CFD方法得到的气动特性数值模拟结果，CRM-WBHS_FSC代表采用流固耦合方法得到的气动特性数值模拟结果，同时图中给出了ETW风洞试验的测力结果。由图可见，采用流固耦合方法，考虑模型的静气动弹性变形影响导致相同迎角下升力系数、阻力系数下降，俯仰力矩系数增加。CRM-WBHS_FSC的模拟结果更接近试验结果。气动特性数值模拟结果与试验结果的差异需要从洞壁干扰、转捩模拟等方面进一步开展工作。

图9 CRM 构型机翼最大变形量收敛历程Fig.9 Convergence process of wing maximal deformation for CRM configuration

(a) CL～α

(b) CD～α

(c) Cm～α

4.2 AGARD445.6机翼颤振数值模拟

采用松耦合方式开展AGARD445.6机翼的颤振特性数值模拟。CFD计算网格规模为2566万，CFD网格的变形采用并行处理方式；结构动力学方程求解取AGARD445.6机翼前四阶模态；CFD计算网格、各阶模态的固有频率如图11所示。图12给出了AGARD445.6机翼颤振计算得到的广义位移响应曲线，来流状态为Ma=0.96，q/qe=0.8650时，计算得到的颤振速度系数为Vf=0.2817，其中qe为试验动压。

图11 AGARD445.6机翼CFD计算网格与结构模态Fig.11 CFD grid and structural modal for AGARD445.6 wing

图12 AGARD445.6机翼广义位移随时间变化历程Fig.12 Time history of generalized displacement for AGARD445.6 wing

图13为AGARD445.6机翼颤振计算得到的各马赫数临界颤振速度系数与试验的比较，马赫数范围0.499～1.072，同时给出了320万网格规模的数值模拟结果。采用320万网格规模进行颤振模拟，采用720个FT1500处理器核，迭代5000个时间步需要65.9 h；采用2566万网格规模进行颤振模拟，采用1600个FT1500处理器核，迭代5000个时间步需要255.6 h。由图13可见，计算结果准确预测了跨声速颤振的“凹坑”现象；在亚跨声速范围内，计算结果与试验结果吻合良好；在超声速区域，计算结果明显高于试验结果，超声速区域计算结果与其他文献的计算结果类似[22-23]。这个具体算例的数值结果一方面验证了本文超大规模颤振问题的数值模拟能力；另一方面，网格规模的增加并没有使得计算结果与试验结果更加吻合。计算结果与试验结果的差异需要从风洞试验以及数值模拟方法两个方面进一步开展工作。

图13 AGARD445.6机翼颤振特性数值模拟Fig.13 Flutter numerical simulation for AGARD445.6 wing

5 结 论

建立了静气动弹性、颤振模拟的一般流程并开发了超大规模气动弹性数值模拟软件的初级版本，实现了百亿量级超大规模静态结构网格生成和十亿量级结构网格自动变形，开展了动网格模块、CFD求解模块、颤振模拟的并行效率和并行一致性测试，采用CRM-WBH构型和AGARD445.6机翼模型进一步开展了气动弹性数值模拟软件的确认工作。

下一步工作主要包括：采用MPI+OpenMP并行策略，进一步提高气动弹性数值模拟的并行效率；开展时域非定常气动弹性计算方法研究，提高非定常气动弹性问题数值模拟精度；利用高保真度静气动弹性、颤振风洞试验数据，进一步开展气动弹性模拟软件的确认工作。

致谢:感谢参加本课题的张书俊、李伟、杨小川等三位同志在超大规模网格生成、结构有限元分析和计算结果处理等方面的工作。