基于多核并行计算的无人机CFD数值模拟

姜悦宁，贾宏光，厉　明

JIANG Yuening1，2，3,JIAHongguang1，3,LI Ming1，3

1.中国科学院 长春光学精密机械与物理研究所，长春 130033

2.中国科学院大学，北京 100039

3.长光卫星技术有限公司，长春 130033

1.Changchun Institute of Optics,Fine Mechanics and Physics,ChineseAcademy of Sciences,Changchun 130033,China

2.University of ChineseAcademy of Sciences,Beijing 100039,China

3.Chang Guang Satellite Technology Co.,Ltd.,Changchun 130033,China

1　引言

无人机以其结构轻便、可操纵性强等独特的技术优势，在航空测绘、摄影和森林植保等领域得到越来越广泛的应用[1]。气动布局设计是无人机设计过程中的关键，通过对全机气动特性进行分析，可以用于进一步分析飞机的飞行性能，为控制系统设计提供重要依据[2-4]。计算流体力学（Computational Fluid Dynamics，CFD）是指通过计算机利用数值方法求解流体动力学方程，获得流动规律并解决流动问题。利用CFD技术可以模拟流体的真实状况，快速获得气动特性数据[5-6]。

航空工程领域是CFD的最大推动者之一，随着计算精度要求的提高，与此同时需要更大规模的计算网格，传统串行的计算方法已经难以满足需求，必须采用并行计算的方法解决大规模流场数值计算耗时过多的问题。多核并行计算是指同时对多个任务、多条指令或多个数据项进行处理[7]。文献[8]对比了多核并行集群系统与单核集群的计算时间，验证了多核集群的高效性和可扩展性；文献[9]基于以太网连接的工作站机群，采用FLUENT并行计算方法对某型号飞机的非结构网格进行了数值求解，并讨论了计算节点数对计算时间、加速比和并行效率的影响；文献[10]采用CFD并行计算的方法对乘波飞行器前体模型进行了计算。然而基于多节点的计算集群受节点间通信的限制，加速比难以达到理想值[11]。

随着硬件技术的快速发展，单机多核CPU应运而生，它将几个甚至几十个CPU集成到单个芯片上，每个CPU核作为一个单独的处理器，从而进一步提高并行计算的速度。

为了提高飞机总体设计中气动布局设计的效率，需要获得一种高效的数值模拟方法，在提高数值计算精确度的同时缩短计算时间。本文首先采用串行的计算方法分别对基于Euler方程和N-S（Navier-Stokes）方程的求解器进行验证，将求解结果与ONERA M6机翼[12]实验结果进行对比，以验证主控方程对计算精度的影响。接着采用FLUENT并行的计算方法，讨论了线程数对M6机翼流场数值计算时间和加速比的影响，寻找一种基于FLUENT的高效并行求解方法。最后得出一种准确、高效的计算方案，并针对某V型尾翼、上单翼布局无人机的整机气动外形进行了计算。

2CFD并行计算技术

CFD技术与理论分析、实验研究共同组成空气动力研究的三种主要手段，相较于其他两种方法具有研制周期短、损耗小的特点，同时不需考虑洞壁干扰和支架干扰，可以在计算机上实现各种条件下的流场计算。随着计算机的飞速发展，使用计算机对工程现象进行数值模拟的分析手段成为可能，并且计算精度和效率都大大依赖于计算机的硬件设备。

CFD任务能够分解为可并行计算的多个子任务，在不同处理器上同时处理，这就为并行计算提供了可行性。采用并行计算一方面能够解决串行计算时间长的问题[13-15]，另一方面可以在相同计算时间内扩大求解规模，提高计算精度[16-18]。

2.1　CFD数值模拟的控制方程及离散方法

Navier-Stokes方程组的求解是CFD的核心，该方程组是典型连续介质流体力学的控制方程。Navier-Stokes方程组包括质量守恒、动量守恒和能量守恒方程，直角坐标系下，积分形式的三维守恒型N-S方程形式如下：

式中Q=[ρ,ρu,ρv,ρw,ρe]T表示守恒向量，n为边界外法向量，∂V为某一固定区域V的边界。

分别表示无粘项和粘性项，ρ、e分别为密度和单位质量气动总能量，(u,v,w)为直角坐标系下的速度分量。其中粘性应力项：

无粘项空间离散采用迎风通量差分分裂格式，对流项和压力项可以通过一个网格单元的通量平衡来表示：

其中UL、UR是交界面的原始参数。

粘性项采用二阶差分格式进行离散：

时间离散采用近似因子分解法（Approximate Factorization Method），该方法是一种基于多维非线性控制方程的隐式方法，从而提高总体计算效率。

2.2　多核并行技术

不同于单核CPU多线程的交错执行，多核CPU的多线程在物理上是并行执行的。要充分发挥多核CPU的硬件性能，必须采用多线程执行，使得每个CPU核在同一时刻都有线程在执行。

ANSYS FLUENT CFD（以下简称FLUENT）中的并行处理包括FLUENT、主机进程和一套节点进程之间的交互作用，通过CORTEX模块与主机进程及各个节点进程进行信息传递。串行求解器各模块间逻辑关系如图1所示。

图1　串行求解器各模块间的逻辑关系

FLUENT的并行计算可在大型并行计算机上运行，也可在一个多CPU的工作站上运行。CORTEX发出的指令通过主机进行解释，由主机进程通过指定计算节点的数据通讯模块传递给其他节点，被指定的计算节点称为“计算节点1”，其他节点在接收主机进程传递的命令后，在各自的数据集合上同时执行相同的程序。各计算节点彼此相连，节点间数据传输、同步和执行全局任务则需要通过数据通讯模块实现，从而实现逻辑意义联通。FLUENT并行求解器各模块间的逻辑关系如图2所示。

图2　并行求解器各模块间的逻辑关系

MPI（Message Passing Interface）是一种信息传递编程标准，为基于信息传递的并行程序提供一个高效统一的编程环境，是目前应用最为广泛的并行编程方式。FLUENT自带的MPI并行机制不仅支持对称多处理共享存储并行计算，也支持网络分布式并行计算，由于多核CPU属于共享存储并行机，因此，采用FLUENT内置MPI并行机制能够有效提高单机多核并行计算的效率。

3　计算应用与分析

3.1　计算条件与并行环境

为了获得高效准确的计算方法，首先对算例机翼进行数值求解，验证求解器的准确性。数值计算采用基于密度的耦合隐式算法，控制方程在计算域的离散方法采用有限体积方法，机体表面满足无滑移边界条件，远场为自由可压缩流场，湍流模型为剪切滑移SST（Shear-Stress Transport）模型。

由于计算工况较多，网格数量较大，普通计算机串行计算将耗费巨大的时间，因此需要并行计算进行实验。实验采用的多核工作站配置两个2.6 GHz的CPU，内存为128 GB，48核数运行处理器。

3.2　串行与并行求解器

本文选取用于验证计算方法的ONERA M6实验机翼作为验证算例，计算条件如表1所示，采用串行的计算方法。

表1　ONERA M6机翼来流属性及计算状态

气动网格为ANSYS ICEM CFD软件生成的结构网格，空间分区及机翼表面网格如图3所示。分别采用N-S方程和Euler方程对流场进行了数值求解，并与实验数据进行了对比。沿机翼展向截面的压力系数如图4所示，可以看出N-S解与实验值吻合很好，这是由于该状态下激波的作用导致流动分离，N-S解由于粘性的计入能够更准确地描述该状态下截面的附体流动。总之，这种三维结构网格划分方法以及流场解算器效果是很好的，可以作为无人机流场数值求解的方法。

图4　ONERA M6亚音速状况下的压强分布（z l=0.2）

加速比是衡量并行计算效果的核心指标，加速比通常是指在相同规模下串行执行时间与并行执行时间的比值。若Tp表示并行系统下的计算时间，Ts表示串行系统计算时间，那么并行化加速比S可以用下述公式来表示：

分别采用串行和多核并行系统对模型进行计算，设置一组FLUENT并行计算的线程数目，分别为1、6、12、24、36、48，表2、表3列出了1 000个时间推进步的计算时间和计算加速比。由表2可以看出，并行计算可以大大缩短运算时间，随着线程数目的增加，计算效率得到进一步提高，当多线程所需的开销增大到一定程度则会降低运算速度。随着网格数量的增加，多线程并行时加速比普遍略有所降低。由此可见，合理选择并行计算的线程数，可以充分发挥并行计算的优势，并用于大规模的求解问题中。由图5可以看出，不同线程的多核并行计算的压强分布吻合一致，且与实验值吻合较好。

表2　多核并行计算与串行计算执行时间对比

表3　多线程数并行计算加速比

3.3　数值计算与分析

网格生成与数值模拟的结果密切相关，网格质量的好坏直接关系到数值模拟结果的正确性。本次实验前处理模型仍采用ANSYS ICEM CFD软件生成计算域内结构网格。由于计算模型的对称性，首先对半区几何模型进行分区结构网格划分，通过网格镜像的功能生成全机网格。

数值计算内容为某测绘无人机无动力状态下的气动力参数，采用上单翼和V型尾翼布局形式。计算条件为：Ma=0.065，Re=4.87×105，选择三维可压缩N-S粘性方程作为流场主控方程，采用FLUENT软件进行多核并行数值计算，通过设定线程数目，提高计算效率。求解时利用软件自带的脚本记录功能，将不同飞行工况计算的步骤记录下来，通过读入脚本实现不同工况自动计算和保存，并在计算结束后自动切换下一工况。这种方法避免了计算资源的浪费，同时减少因切换工况的时间花费。图6为整机升力系数和阻力系数随迎角的变化曲线；图7为机体表面的等值线分布和表面压力系数分布云图。

4　结论

（1）本文提出了一种用于求解无人机整机的精确CFD数值模拟方法。分别采用Euler法和N-S法对算例机翼进行解算，验证了以N-S方程作为流场主控方程的准确性，并采用该解算器对某型号民用无人机流场的气动特性进行了数值模拟。

图5　ONERA M6机翼不同线程数各截面的压强分布

图6　整机气动特性曲线

图7　机体表面等值线分布和压强分布云图

（2）提出了一种用于求解无人机整机的单机并行计算方法。本文基于多核系统对算例机翼的求解执行时间进行了探讨，验证了并行求解的准确性和高效性，通过设置不同线程数目执行求解程序，使CPU不同核心同时进行计算，达到加速的目的。随着线程数的增加，加速比的提升幅度先大后小，达到最大值后会有所减小。将这种思路用于大规模的无人机整机数值模拟问题中，使得计算效率得到大幅度提高。

（3）由于多核处理器多个CPU核的共享存储模式，相对于集群计算平台，节点间的通信时间几乎可以忽略，可同时兼顾计算机硬件结构和发挥CPU最大计算能力，具有很大工程应用潜力。进一步的研究内容包括：流场网格分区数与计算节点数的关系，即如何对三维流场进行分区，并选择合适的计算节点，避免硬件资源的浪费。

参考文献：

[1]崔秀敏，王维军，方振平.小型无人机发展现状及其相关问题分析[J].飞行力学，2005，23（1）：14-18.

[2]祝小平，向锦武，张才文，等.无人机设计手册[M].北京：国防工业出版社，2007.

[3]段镇，高九州，贾宏光，等.无人机滑跑线性化建模与增益调节纠偏控制[J].光学精密工程，2014，22（6）：1507-1516.

[4]李迪，陈向坚，续志军.增益自适应滑模控制器在微型飞行器姿态控制中的应用[J].光学精密工程，2013，21（5）：1183-1191.

[5]关键.低速无人机动态气动特性数值模拟及布局研究[D].长沙：国防科技大学，2013.

[6]Baker A M，Ying S X.A fixed time performance evaluation of parallel CFD applications[J].Super Computing，2002，95：18-23.

[7]周伟明.多核计算与程序设计[M].武汉：华中科技大学出版社，2009.

[8]陈坚祯，阳平，李斌，等，多核并行计算下的流量传感器流场模拟研究[J].衡阳师范学院学报，2011，32（6）：82-84.

[9]邵波，毛国勇，张武.基于Fluent的全机数值模拟及并行计算[J].计算机工程与设计，2006，27（17）：3178-3180.

[10]潘沙，李桦，夏智勋.高性能并行计算在航空航天CFD数值模拟中的应用[J].计算机工程与科学，2012，34（8）：191-198.

[11]熊玮，夏文龙，余晓鸿，等.多核并行计算技术在电力系统短路计算中的应用[J].电力系统自动化，2011，35（8）：49-52.

[12]Schmitt V，Charpin F.Pressure distributions on the ONERAM6 wing at transonic mach numbers，AGARD-AR-138[R].1979.

[13]冯定华，潘沙，丁国昊，等.CFD并行计算方法在高超声速数值模拟中的应用[C]//第三届高超声速科技学术会议，2010.

[14]吴颀.GPU并行计算及其在飞行器设计中的应用[D].北京：北京理工大学，2015.

[15]唐逸豪，高振勋，蒋崇文，等.基于LPT近似算法的CFD并行计算网格分配算法[J].工程力学，2015，32（5）：243-249.

[16]Anderson J D.Computational fluid dynamics[M].[S.l.]：Mc Graw Hill Education，2007.

[17]刘巍，张理论，邓小刚.计算空气动力学并行编程基础[M].北京：国防工业出版社，2013.

[18]温建华，朱自强，吴宗成.基于N-S方程串并行计算的机翼优化设计[J].北京航空航天大学学报，2008，34（2）：127-130.