基于Euler方程的返回舱气动外形优化设计方法研究

李治宇，杨彦广，袁先旭，唐志共

（中国空气动力研究与发展中心，四川 绵阳 621000）

0 引 言

气动外形的优化设计是飞行器设计的重要环节。由于飞行器方案论证阶段气动外形优化设计是通过数值优化方法和气动性能计算的结合来实现的，其可靠性、精度和效率是关系到设计成败的重要因素。而精度和效率主要由流场解算和优化算法决定。当前的工程算法求解速度快但精度较低，实用价值有限。在优化设计方面，气动外形的优化设计一般采用已趋于成熟的梯度法等优化方法，但其最大的缺点在于设计结果容易陷入局部最优，并且对目标函数和约束的要求比较高。模拟退火算法等全局优化算法的计算量比较大，应用NS方程求解的时间长。而多目标优化设计经常选用的加权和方法，其权值不容易合理确定，并且不能很好地反应各优化目标之间的关系。

国内外研究者们针对各类飞行器外形进行了大量优化设计工作。印度的Rajesh Kumar Arora等人用牛顿近似的方法求解流场载荷，在法线方向引入压力修正，通过多目标遗传算法（Multi Objective Genetic Algorithm）优化设计了可重复使用的飞船返回舱外形，得到了满足返回舱气动性能并且可重复使用、成本较低的外形［1］。NASA的James L.Brown等人通过优化设计非对称防热底来提高升阻比、降低对流和辐射传热，其中气动力计算部分调用了高保真度工程算法的CBAero软件包，用修正牛顿理论方法计算表面压力分布［2］。Alan Le Moigne和 Wai Sam Wong等人通过求解N－S方程计算翼型的气动性能，并应用序列二次规划（NLPQL）方法进行优化设计。该方法能够快速得到优化解，精度较高，但却容易陷入局部最优的境地［3－4］。唐伟、车竞、夏露等人应用内伏牛顿理论及定常势流格林函数法计算飞行器气动性能，并且通过遗传算法优化设计了新外形的气动布局［5－9］。

本文基于iSIGHT优化软件平台，应用遗传算法和数值优化算法相结合的混合算法，通过应用快速求解Euler方程的CART3D程序，发展了一种快速有效的气动布局优化设计方法，并以CTV返回舱为例验证该方法，将其应用于类神舟飞船返回舱的优化设计中得到可靠结果。一般说来，该方法可用于飞行器气动布局的优化设计。

1 研究方法

1.1 研究方法简介

通过SolidWorks软件作图得到基本气动布局外形并完成参数化建模。在优化设计过程中用到的参数化建模方法适用于SolidWorks和CFD分析软件。本文应用一个VBS脚本文件来完成外形的参数化，该VBS文件通过直接修改SolidWorks文件的特征尺寸来修改外形，并输出parasolid格式的数模文件。在输入错误尺寸数据时，如负值，该文件并不会中断程序，而是在输出原始外形后继续优化过程，优化过程忽略错误数据。所以，在出现错误尺寸时该文件不会导致优化过程中断，也不会影响到优化结果，具有很好的鲁棒性。

应用快速气动分析软件CART3D计算外形的气动性能。该软件通过快速生成非结构笛卡尔网格、求解Euler方程的方法计算高超声速气动力特性。在第一次运算中记录操作命令流文件，之后的气动计算就可以直接通过所生成的命令流文件来运行CART3D程序，期间不需要人工干预操作。并且，iSIGHT可直接调用该文件并修改其中参数。该软件精度高、求解速度快、网格生成高度自动化、可用于机翼、全机、导弹及喷管等复杂系统的气动力计算分析，适用性很好［10］。

优化计算过程选择iSIGHT优化平台，调用合适的优化算法，完成优化设计部分的工作，最终得到优化结果。iSIGHT平台将大量需要人工完成的工作实现自动化处理，并且易于处理复杂问题的优化设计过程，可采用多种方法探索，可缩短设计周期、降低设计成本，提高设计质量和可靠性。iSIGHT提供了多种商业软件及自编程序接口，设计者通过图形用户界面GUI可实现设计问题的过程集成、问题表述、优化方案选择及求解监控等，iSIGHT提供多种优化算法，通过组合或单独使用可基本满足设计者要求。

集成的整个优化过程包括四部分，分别为外形特征尺寸计算（calculation0）、参数化生成外形（SW）、网格生成及气动力计算（CART）、目标函数计算（calculation1）。如图1所示。

图1 优化流程图Fig.1 Flow chart of the optimization

1.2 流场求解方法验证

1.2.1 解算精度验证

CART3D求解器基于可压缩Euler方程，采用有限体积法离散，流动变量位于网格中心，时间推进采用Runge－Kutta法得到稳态解，空间离散为迎风格式，选择使用限制器，格式具有TVD性质，采用多重网格法加速迭代过程。

通过CART3D计算得到飞船返回舱的高超声速气动力系数，并与文献数据对比，验证计算方法的精度。计算来流马赫数M＝6，攻角α＝20°，图2为计算网格和流场压力分布图，表1为计算结果对比，表中数据误差为CART3D计算结果与各参考文献数据对比的相对误差，结果显示该方法计算精度是可以信赖的。其计算质量主要取决于所生成的笛卡尔网格的大小，可以通过调整网格加以改进。图3为CART3D和FLY3D计算的飞船返回舱升阻比和对顶点的俯仰力矩系数随攻角变化的曲线及相应实验结果，可见其计算结果符合很好。

图2 计算网格与流场压力分布图Fig.2 Cartesian grid and flow－fields pressure distribution

表1 计算结果与文献数据对比Table 1 The numerical results vs.the reference

图3 返回舱升阻比和俯仰力矩随攻角变化的曲线Fig.3 Comparison of aerodynamic coefficient CL／CDand Cmvs.angles of attack

1.2.2 解算效率验证

通过与FLY3D程序对比发现，FLY3D网格约190万时，在PC机上运行一次的计算机时约6h；而当CART3D网格量为130万时，在PC机上运行一次的时间小于10min。通过表2对比，机时与网格量有类似正比关系，而计算结果的精度在一定范围内变化不大。这样，我们就可以在保证计算精度的前提下适当减少网格，以缩短机时，使得CART3D的计算效率满足优化设计多轮计算的要求。在优化过程中CART3D实际网格量约在20万左右，在PC机上运行一次的时间约为70～90s，满足优化设计要求。表2中计算条件与表1相同。另外，其网格生成高度自动化，仅需数秒时间，效率很高。

表2 不同网格量的CART3D运行时间与精度对比表Table 2 Comparison of run time and precision vs.grid quantity of CART3D

2 优化设计算例验证

2.1 计算模型

［13］，通过完成欧洲航天局为其载人天地往返运输系统计划中宇航员输送舱（Cosmonaut Transportation Vehicle，以下简称CTV）所提出的气动构型CTV的优化设计来对比选择合适的优化算法。CTV构型为钝球双锥外形，其基本构型特征尺寸及轮廓如图4所示。

图4 CTV返回舱轮廓线Fig.4 Contour line of CTV reentry capsule

2.2 计算目标及约束

为确定CTV的最优气动布局，这里采用以下简化模型。在底部直径D和总体积V不变的情况下，满足配平升阻比CL／CD不小于0.8的约束，设计目标为最大的体积利用率RV。设计变量为球头半径RN、前锥长LF、前锥半锥角θF。可表示为：

设计目标： minU＝1／RV

约束条件：D＝3.5m

V＝3m3

CL／CD＞0.8

设计变量： 0.5m＜RN＜1.0m，2.0m＜LF＜3.0m，10.0°＜θF＜30.0°

其中，体积利用率RV＝6／S3／2，V为体积，S为表面积。设计点来流马赫数M＝6，配平攻角αT＝20°。

2.3 算法选择

本文通过不同算法进行优化对比。选用算法包括数值优化算法：MMFD（Modified Method of Feasible Directions－修正可行方向算法）、NLPQL（序列二次规划算法）、LSGRG（Generalized Reduced Gradient－广义简约下降梯度法）及全局优化算法MIGA（Multi－Island Genetic Algorithm－多岛遗传算法）。

数值优化算法通过梯度搜索，处理不等式和等式约束，快速得到最优解；MIGA通过选择、交叉及变异等遗传操作寻优，把种群分组提高全局搜索能力。

2.4 结果分析

所选用的不同算法所得到的优化解基本达到了优化目的，即在满足约束的条件下，优化目标体积利用率RV达到最大值，并且除算法2因为在局部最优点中止优化外，其余算法得到的结果，都略优于文献数据。对比结果如表3所示，表中长度单位为米，角度单位为度。

由表3中算法1与2对比可见，MMFD算法得到的是局部最优解。对比算法3、4、5发现，最优值接近。由优化数据可以得到，算法3中的NLPQL实际上没有起到优化作用，即该算法的优化部分没有得到更优解。算法4的MMFD及算法5的LSGRG的运行对结果进行了更进一步的局部优化。可以说明的是MIGA在不同初始点开始计算，得到了近似的全局最优点。而从算法6中可以看到，将MIGA算法的进化代数由10改为12后，得到了“更优解”。所得最优解的气动力系数如表4所示，与文献［14］中配平攻角20°时最大升阻比为0.8左右的计算结论及FLY3D精细计算的结果相吻合。

表3 不同优化算法的优化结果对比Table 3 Comparison between different optimization techniques

表4 算法6得到最优解的气动力系数Table 4 Aerodynamic coefficient of the optimization result for algorithms 6

通过以上分析，可以确定适合于气动布局外形优化的优化算法。通过多岛遗传算法MIGA全局优化，得到全局最优解后再通过数值优化算法进行优化计算，得到局部更优解，而算法MIGA的参数设置将会影响到结果的优化程度。这样的混合优化算法在避免陷入局部最优解的同时，数值优化算法的优化计算对全局优化算法得到的最优解进行“微调”，得到“更优解”。

3 类神舟飞船返回舱优化设计

3.1 单目标优化设计

类神舟飞船返回舱外形的主要特征参数包括最大截面直径dm、球冠半径Rn、大头肩部半径Rc、倒锥角θ、长度L、尾部圆弧半径RB及尾部圆台的底面半径和高度。如图5所示。

图5 类神舟飞船返回舱轮廓线Fig.5 Contour line of Shenzhou Spaceship reentry capsule

将各特征参数，除倒锥角外，用dm无量纲化。球冠钝度大头肩部圆弧曲率半径Rc／dm、尾部圆弧曲率半径长细比L／dm。优化要求在体积不小于16m3、体积利用率不小于0.85时，通过变化最大截面直径、球冠钝度、大头肩部曲率半径及倒锥角四个变量得到升阻比最大的优化外形，计算来流M＝13.2，αT＝20.0°，优化计算模型表示为：

目标： maxCL／CD

约束：V≥16m3

Rv≥0.85

变量： 2.5m≤dm≤4.0m

选用MIGA和MMFD相结合的优化算法，其中的多岛遗传算法设置岛数为10，每岛10个个体，进化12代，染色体长度20，交叉概率0.9，变异概率0.01，迁徙率0.4，迁徙间隙3代。

3.2 双目标优化设计

双目标优化设计与以上单目标优化设计相比，约束条件及变量均相同，不同之处在于设计目标要求在升阻比最大的同时峰值热流最小化，其中峰值热流的计算采用Fay－Riddell公式。优化计算模型可表示为：

表5 单目标优化设计结果Table 5 Optimization result for single objective

目标： maxCL／CD，minqw

约束：V≥16m3

Rv≥0.85

变量： 2.5m≤dm≤4.0m

选用邻域培植遗传算法NCGA，种群大小20，进化150代，最后输出pareto解集。两个目标的比例因子分别取为1.0和1000.0，权系数均为1.0，则目标函数objective＝qw／1000.0－CL／CD／1.0。相应优化结果及优化外形如表6、图6所示。

表6 双目标优化设计结果Table 6 Optimization result for double objectives

图6 类神舟飞船返回舱优化pareto最优解集及最优化外形图Fig.6 Pareto plot and the optimization contour of Shenzhou Spaceship reentry capsule

4 小 结

（1）形成了基于iSIGHT软件平台的集成优化方法，用于返回舱优化设计能够得到理想结果。其优化精度高、速度快，有较好的工程应用前景。

（2）所采用的流场解算方法适用于多种复杂系统的气动力计算，适用性强，计算精度和效率较高。

（3）选择MIGA与数值优化算法相结合的混合遗传算法，可以满足单目标多变量优化设计。NCGA同样可以满足多目标优化设计，影响优化结果的主要因素是遗传算法选择、交叉及变异等相关参数的选择。下一步的研究工作将集中在优化算法的选择应用及相应优化参数对优化结果的影响方面。

参考文献:

［1］RAJESH KUMAR ARORA.Aerodynamic shape optimization of a Re－entry capsule［R］.AIAA 2003－5394.

［2］JAMES L BROWN，JOSEPH A GARCIA.An asymmetric capsule vehicle geometry study for CEV［R］.AIAA 2007－604.

［3］ALAN LE MOIGNE.A discrete Navier－Stokes adjoint method for aerodynamic optimisation of Blended Wing－Body configurationsp［D］.［PHD Thesis］.Cranfield U－niversity，2002.

［4］WAI SAM WONG.Mechanisms and optimisations of 3Dshock control bumps［D］.［PHD Thesis］.Sheffield University，2006.

［5］唐伟，桂业伟，张勇，等.未来升力再入飞船返回舱的气动布局设计［C］.中国宇航学会首届学术年会.北海，2005.

［6］唐伟，高晓成，李为吉，等.双椭圆截面再入飞行器的气动计算及布局优化设计［J］.空气动力学学报，2004，22（2）：171－174.

［7］唐伟，桂业伟，方方.新型升力再入飞船返回舱气动外形选型研究［J］.宇航学报，2008，29（1）：84－88.

［8］车竞，唐硕，何开锋，等.高超声速飞行器气动布局总体性能优化设计研究［J］.空气动力学学报，2009，27（2）：214－219.

［9］夏露，高正红，李天.飞行器外形多目标多学科综合优化设计方法研究［J］.空气动力学学报，2003，23（1）：275－280.

［10］JAMES S GREATHOUSE，BENJAMIN S KIRK，RANDOLPH P LILLARD.Crew exploration vehicle（CEV）crew module shape selection analysis and CEV aeroscience PROJECT overview［R］.AIAA 2007－603.

［11］袁先旭.非定常流动数值模拟及飞行器动态特性分析研究［J］.［博士学位论文］.中国空气动力研究与发展中心，2002.

［12］纪楚群，周伟江.钟形返回舱空气动力特性数值模拟［J］.航天返回与遥感，2001，22（1）：32－37.

［13］安伟刚.多目标优化方法研究及其工程应用［D］.［博士学位论文］.西北工业大学，2005.

［14］唐伟，张勇，李为吉，等.欧洲返回舱CTV气动特性综述［J］.空气动力学学报，2003，21（3）：362－366.