引言

随着航空燃油价格的不断上升，改善燃油经济性是现代飞机设计追求的目标，采用展弦比更大的薄机翼有利于进一步提高飞机的气动特性，然而机翼柔性的增大将会造成更为明显的气动弹性变形，对飞机气动特性造成显著的影响^[1]。但是传统的设计方法在进行气动外形设计时，并没有很好的考虑静气动弹性的影响，为了维持有利的巡航外形，往往采用增大机翼刚度的方法来保证飞机具有足够的气动弹性稳定性，结构强度以及避免不同飞行状态和机动操纵时机翼发生不利的扭转变形。而这种方式将会显著增加机翼的结构重量。因此随着复合材料等新材料的大量使用，以及飞机柔性的增大，需要在初步设计阶段就考虑柔性机翼的气动/结构耦合问题。

对于柔性较大的飞机，不同的飞行状态对应不同的载荷分布，进而影响飞机的静气动弹性变形, 而静气动弹性变形产生的扭转角分布变化以及机翼明显的上弯扰度变形所形成的弧形机翼^[2-3]对飞机的气动特性会带来明显的影响。柔性机翼气动/结构一体化设计充分利用了机翼的柔性特点，结合气动弹性分析技术，通过合理的气动、结构耦合设计使得在多种飞行状态下机翼都能够获得期望的气动弹性变形，以提高飞机的综合性能。

国内外学者在气动/结构一体化设计方面都做了大量的研究工作。董波等^[4]将全速势方程CFD数值评估方法与基于工程梁理论的机翼结构设计方法相结合，进行了气动/结构一体化设计，但是没有考虑静气动弹性的影响。余雄庆^[5]等利用代理模型预测气动力系数和翼尖位移、扭转变形，以机翼剖面翼型扭转角和结构尺寸为设计变量进行考虑静气动弹性影响的气动/结构一体化设计的研究，并对设计结果进行了简单的描述。代理模型的采用在一定程度上提高了计算效率。陈海昕^[6]等利用智能算法采用“三步”多目标气动/结构优化设计策略对机翼平面形状、机翼扭转角、翼型剖面外形和结构尺寸参数进行了气动/结构多学科优化设计研究。“三步”优化策略在一定程度上将气动，结构进行解耦，减少了具有大规模设计变量的复杂气动/结构多学科优化问题的计算量。但是，这种优化策略也在一定程度上削弱了气动、结构之间的耦合作用。左英桃^[7]利用离散伴随方法将基于欧拉方程的气动力模型和有限元结构求解器相结合，进行了以剖面翼型和扭转角为设计变量的，考虑了静气动弹性影响的气动优化设计研究。伴随方法的使用显著地提高了优化设计的效率。经过优化设计，全机升阻比提高了4.06%，但是整个优化设计过程并没考虑结构参数的变化。

在国外，以美国Martins团队为代表的课题组在以伴随方法为基础的气动/结构一体化设计方法方面的研究较为成熟。Martins^[8-9]利用伴随方法建立了基于RANS方程和TACS (the Analysis of Commposite Structures)结构求解器的针对复杂构型大设计变量的气动/结构优化设计框架，并成功地应用到了大展弦比跨声速客机机翼的设计研究中。

虽然目前国内外在气动/结构一体化设计方面有大量的研究可寻。但是，大多数研究工作都是围绕着以提高气动/结构一体化设计方法的设计效率、计算精度和计算复杂度开展。少有论文从应用角度出发，将气动/结构一体化设计与单学科的优化设计结果进行对比研究，较详细地分析考虑静气动弹性影响之后，对机翼气动、结构设计带来的影响。

因此，本文从应用角度出发，针对现代大型跨声速民用客机开展气动/结构一体化设计问题的研究。在设计方法上，为了兼顾计算效率和计算复杂度之间的矛盾，采用基于全速势方程加附面层修正的CFD数值评估方法，板壳有限元模型和松耦合静气动弹性分析流程的气动/结构一体化设计方法。针对0.85马赫数跨声速客机的机翼设计，利用建立的优化设计方法进行了分别以机翼扭转角分布和剖面翼型为设计变量的单点优化设计研究。并将气动/结构一体化设计结果与单学科的气动设计结果进行对比，初步研究了考虑静气动弹性影响后，对机翼气动载荷、翼型厚度和弯度分布，以及结构尺寸厚度分布设计等方面带来的影响。

1 全速势方程加附面层修正的气动数值评估方法

随着CFD技术在飞机气动外形设计中的广泛应用，数值评估方法的可靠性、准确性及计算速度决定了设计阶段的效率。研究^[10-12]表明飞机在巡航状态没有强激波及大分离区存在时，是可以采用全速势方程进行流场的快速评估的，而对于更为复杂的流动现象，需要更高精度的数值模拟方法。本文采用的全速势流场评估方法通过在边界层理论范围内建立有黏-无黏耦合迭代进行流场求解。首先由全速势方程求解全流场，得到流场的无黏解，再以无黏解为边界条件，由可压缩附面层方程计算附面层流动，得到附面层流动，以及附面层内各流动参数。然后再次计算全速势方程，如此进行有黏无黏迭代求解直到收敛。

跨声速无黏外流场采用守恒形式的全速势方程，利用有限差分法进行求解，全速势控制方程如下所示：

$ {\left( {\frac{{\rho U}}{J}} \right)_\xi } + {\left( {\frac{{\rho V}}{J}} \right)_\eta } + {\left( {\frac{{\rho W}}{J}} \right)_\zeta } = 0 $

(1)

$ \left( \begin{array}{l} U\\ V\\ W \end{array} \right) = \left( \begin{array}{l} {\phi _\xi }\\ {\phi _\eta }\\ {\phi _\zeta } \end{array} \right) $

(2)

$ \left\{ \begin{array}{l} {\phi _\xi } = {({\phi _{i + 1}} - {\phi _i})_{j,k}}\\ {\phi _\eta } = \frac{1}{4}{({\phi _{i + 1,j + 1}} - {\phi _{i + 1,j - 1}} + {\phi _{i,j + 1}} - {\phi _{i,j - 1}})_k}\\ {\phi _\zeta } = \frac{1}{4}{({\phi _{i + 1,k + 1}} - {\phi _{i + 1,k - 1}} + {\phi _{i,k + 1}} - {\phi _{i,k - 1}})_j} \end{array} \right. $

(3)

通过可压缩附面层方程和有限差分格式进行黏性附面层求解。

通过将全速势方程加附面层修正的数值方法应用于标模(Common Research Model, CRM)^[13]在设计点(Ma=0.85, C_l=0.5, Re=1×10⁶)的计算，并与实验数据和RANS方程(采用SST湍流模型，Roe格式，网格量800万)计算结果进行对比，验证所采用的数值方法的可靠性。CRM见图 1。

由图 2可以看出，在沿展向的6个剖面位置处，全速势方程加附面层修正的计算结果与RANS的计算结果总体比较吻合，并与实验数据趋势一致，尤其头部峰值和下翼面的压力分布重合度比较高，两种数值计算方法得到的压力恢复梯度、激波位置和激波强度略有差异，但是整个机翼的压力云图(图 3所示)形态一致，在工程误差范围之内。因此，在初始机翼气动力设计阶段，在没有强激波和大分离流动时，完全可以利用全速势方程加附面层修正的方法代替精度较高的RANS方法，作为马赫数0.85的高马赫数超临界机翼的气动评估手段，以显著提高设计效。

2 静气动弹性计算方法及算例验证 2.1 静气动弹性分析流程

求解气动弹性问题的耦合方法通常可以分为紧耦合和松耦合两种。松耦合方法以模块化的形式进行气动弹性计算，克服了紧耦合方法面临的计算求解困难的问题^[14]。松耦合方法气动力模型和结构模型不断耦合，交换数据，当两者达到收敛时，得到一个最终的静气动弹性分析结果。本文采用C2型基函数的RBF插值方法^[15-16]作为气动与结构求解器之间的数据交换技术，全速势方程加附面层修正的方法作为气动力数值评估方法。整个静气动弹性分析流程如图 4所示。

2.2 DLR F6翼身组合体验证算例

选取DLR-F6翼身组合体构型进行静气动弹性分析程序的验证。计算状态为Ma=0.75, Re=3×10⁶, 定升力系数C_L=0.5。机翼有限元模型采用实体模拟机翼，并挖空其中的翼盒部分，通过MSC PATRAN进行有限元建模。材料弹性模量E=2.05×10¹¹, 泊松比υ=0.3。

图 5和图 6分别为本文所建立的静气动弹性分析方法与用Workbench商用软件计算得到的沿展向Z方向后缘挠度变形、扭转变形与实验数据^[17]的对比图。可以看出两种不同计算方法得到的Z向后缘挠度变形和扭转角变形基本贴合，仅在翼梢附近具有一定明显差别。与实验数据相比，两种不同的计算方法得到的变形结果整体趋势与实验数据较为吻合，但仍存在一定误差。

对于静气动弹性，结构变形与时间无关，因此静气动弹性分析取决于结构刚度分布的准确性、定常气动载荷计算的精度和数值交换的准确性。图 7为分别采用全速势方程加附面层修正的数值方法与RANS方程计算(采用SST湍流模型，Roe格式，网格量800万)F6标模得到的沿展向4个站位处的剖面压力分布结果以及与实验数据对比图，从图中可以看出两种CFD数值评估方法计算得到的压力分布除在翼梢附近及头部峰值外基本重合，并与实验数据贴合相对较好，说明所采用的全速势方程加附面层修正的数值评估方法能够比较准确的给出机翼的载荷分布和剖面压力分布。由于无法得到完全符合实验模型的结构有限元模型，因此两种不同计算方法得到的计算结果与实验数据的误差主要来源于结构有限元模型与实验模型之间的差异。计算结果表明所建立的静气动弹性分析方法是可靠的。

3 多学科优化设计框架

采用优化设计流程图表示方法^[18]，进行多学科优化设计框架的描述，如图 8所示。其中灰色粗线为数据流，黑色细线为进程流，X_a(0), X_s(0) 分别为气动和结构设计变量的初始值，X_a^*, X_s^*分别为整个系统得到的最终的气动和结构设计变量。对于考虑静气动弹性影响的气动/结构一体化设计，结构优化设计借助MD NASTRAN，在满足强度和刚度约束的条件下以结构重量最轻为目标, 然后进行静气动弹性的评估，再从静气动弹性计算结果中提取载荷C_p和变形约束进行结构参数的修正。由于结构设计模块的计算量远小于气动模块，因此将结构设计嵌套在气动优化设计中，并最终将静气动弹性分析得到的综合性能反馈给整个系统的优化算法进行全局寻优。以此将结构尺寸设计变量和气动设计变量进行解耦，虽然设计结果的全局性受到了一定的限制，但是大大降低了整个系统的计算量。

整个系统优化采用改进的微分进化算法。1996年微分进化算法参加了首届IEEE进化算法大赛，在所有参赛的进化算法中，DE被证明为最优的进化算法，随后该算法在各个领域都得到了广泛的应用^[19]。DE算法的基本操作包括变异、交叉和选择，标准微分进化算法的进本操作计算如下：

变异算子：

$ \mathit{\boldsymbol{V}}_{_i}^{^{(n + 1)}} = \mathit{\boldsymbol{X}}_{_{r1}}^{^{(n)}} + F\cdot(\mathit{\boldsymbol{X}}_{_{r2}}^{^{(n)}} - \mathit{\boldsymbol{X}}_{_{r3}}^{^{(n)}}) $

(4)

交叉算子：

$ \left\{ \begin{array}{l} u_{_{i,j}}^{^{(n)}} = v_{_{i,j}}^{^{(n)}},{\rm{if}}(ran{d_j}\left[ {0,1} \right] \le CR.\;{\rm{or}}.j = {j_{rand}})\\ x_{_{i,j}}^{^{(n)}},{\rm{otherwise}} \end{array} \right. $

(5)

选择算子：

$ X_{_i}^{^{(n + 1)}} = \left\{ \begin{array}{l} \mathit{\boldsymbol{U}}_{_i}^{^{(n + 1)}},{\rm{if}}(f(U_{_i}^{^{(n + 1)}}) \ge f(X_{_i}^{^{(n)}}))\\ \mathit{\boldsymbol{X}}_{_i}^{^{(n)}},{\rm{otherwise}} \end{array} \right. $

(6)

其中，第n代种群中的第i个个体可表示为X_i⁽ⁿ⁾=(x_{i, 1}⁽ⁿ⁾, x_{i, 2}⁽ⁿ⁾, …, x_{i, D}⁽ⁿ⁾)，D为个体变量的维数；变异操作后得到的中间个体记为V_i⁽ⁿ⁺¹⁾，其中r1, r2, r3∈{1, 2, …, D}，为互不相等的实数；F是一个常数，用来控制变异的比例，通常取F∈(0, 1)；U_i⁽ⁿ⁾=(u_{i, 1}⁽ⁿ⁾, u_{i, 2}⁽ⁿ⁾, …，u_{i, D}⁽ⁿ⁾)为交叉后产生的个体；f(X_i⁽ⁿ⁾)为X_i⁽ⁿ⁾的适应值。

通过将交叉算子中交叉对象按适应值排序，并将次好的个体和差的个体分别与最好的个体进行交叉，同时求取两个交叉算子所得的平均值作为最终交叉结果对基本微分进化算法进行改进，在保证算法全局性的情况下，提高算法的收敛速度。改进后的变异算子为：

$ \begin{array}{l} \mathit{\boldsymbol{V}}_{_i}^{^{(n + 1)}} = [\mathit{\boldsymbol{X}}_{_{r1}}^{^{(n)}} + {F_1}(\mathit{\boldsymbol{X}}_{_{r1}}^{^{(n)}} - \mathit{\boldsymbol{X}}_{_{r2}}^{^{(n)}})]/2 + \\ \;\;\;\;\;\;\;\;\;\;[\mathit{\boldsymbol{X}}_{_{r1}}^{^{(n)}} + {F_2}(\mathit{\boldsymbol{X}}_{_{r1}}^{^{(n)}} - \mathit{\boldsymbol{X}}_{_{r3}}^{^{(n)}})]/2 \end{array} $

(7)

整个系统级的优化以航程最大为目标，为了保证气动学科计算时升力系数不变，将机翼重量的减少量增加为燃油重量，如公式(8) 所示。最后，整个系统以计算得到的航程为依据，将气动和结构设计结果进行综合评估，引导整个系统的优化方向。

$ Range = \frac{{\eta KMa}}{{g{q_e}}}\cdot{a_H}{\rm{ln}}(\frac{{{W_0}}}{{{W_1} + W}}) $

(8)

其中，K为巡航升阻比，q_e为发动机耗油率，W₀为巡航起始阶段飞机重量，W₁为零油重量减机翼结构重量，W为机翼结构重量，表达式为:

$ W = c\cdot{W_1} + {W_2} $

(9)

其中，W₁为直接通过有限元模型计算得到的结构重量；c为放大系数^[20]，用于考虑由于有限元模型的简化而忽略掉的其它部件(不包括增升装置)以及部件之间的链接构件等的重量；W₂为忽略掉的前后缘增升装置及其滑轨等相关部件的重量，因为文中平面形状以及增升装置的布置形式不变，因此为了简化处理，可近似认为这部分重量不发生变化。

4 优化算例 4.1 考虑静气动弹性影响的机翼扭转角优化设计

选取现代大型宽体客机(起飞总重为218000kg)翼身组合体进行考虑静气动弹性影响的机翼气动/结构一体化设计。其中气动设计仅仅考虑扭转角变化，从翼根(机翼安装角不变)到翼梢选取9个控制剖面，结构设计仅仅进行承力结构的尺寸优化。结构有限元模型只考虑了组成承力翼盒的蒙皮、前后梁腹板和翼肋，长桁和梁缘条厚度采用文献[21]和文献[22]中的方法将其“打扁”计入蒙皮厚度，以减小建模复杂度。机翼结构有限元模型如图 9所示，材料为铝材。每两个翼肋作为一个设计变量，前后梁以及上下蒙皮被翼肋划分为多块，其中梁每两块为一个设计变量。上、下蒙皮每一块为一个设计变量。总共有325个设计变量，其中气动设计变量9个，结构设计变量316个，如表 1所示。气动计算状态为Ma=0.85，C_L=0.515，结构设计变量的取值范围为，后梁厚度范围为10~50mm, 前梁厚度范围5~40mm，上、下蒙皮厚度变化范围5~50mm，翼肋厚度变化范围5~30mm。优化目标为航程最大，结构约束包括应力屈服极限(考虑2.5倍的过载系数)约束以及翼梢弯曲和扭转变形约束。为了对比研究考虑静气动弹性影响的多学科优化设计对设计结果产生的影响，进行了以扭转角为设计变量的单学科的气动优化设计，并将初始构型(Ori_aeroelasicity)、考虑静气动弹性影响的多学科设计结果(Opt_aeroelasicity)与单学科气动设计优化结果(Opt_aerodynamic)利用RANS方程对设计结果进行校核。

经过多学科优化，Opt_aeroelascity构型相比于初始构型而言航程提高了5.63%，如表 2所示。其中升阻比提高了6.53%，机翼重量增加了2.56%。初始构型静气动弹性变形后翼梢挠度为1.34m，优化设计结果的翼梢挠度为1.61m(图 10)。从初始构型与优化结果的表面压力云图(图 11)对比可以看出，优化设计结果相比于初始构型载荷适当外移，显著的改善了机翼的气动特性，同时载荷的外移使得机翼结构重量有所增加。由于建立的有限元模型将梁缘条以及长桁的厚度计入蒙皮中，因此弯矩主要由上下蒙皮承受，扭矩主要由梁腹板以及翼盒承受，因此优化结果显示蒙皮厚度普遍较厚。从图 12至图 13可以看出，优化后外翼段载荷较大使得中段翼承受更大的弯矩，因此中段翼蒙皮厚度增加明显，同时载荷的外移使得翼根弯矩增大，因此翼根处蒙皮厚度也有所增大。由于采用的有限元模型没有考虑开口以及由起落架、发动机引起的集中载荷的影响，同时气动仅仅考虑了巡航点的气动性能，因此优化结果与实际机翼结构存在差异。

图 14为单学科气动优化设计(蓝色虚线)与多学科设计(绿色实线)以及初始构型(红色实线)的机翼升力系数和载荷分布对比图。单学科气动优化设计结果载荷分布更加贴近椭圆形载荷分布，因此具有更好的升阻特性。而考虑静气动弹性影响的设计结果载荷分布介于初始构型与单学科气动设计之间，为典型的三角形分布。

图 15为5个典型剖面的压力分布对比图，图 16为扭转角分布对比。可以看出，单学科气动优化设计为了获得较高的升阻特性，减小了外翼段的负扭转，且整个扭转角分布成波浪形。而气动/结构一体化设计结果一方面适当的减小了外翼段负扭转，同时有效的控制了沿展向的扭转角变化，避免出现明显的连续波浪形，除了展向11m位置处略有起伏外，总体上分布比较光滑。可见综合性能最好的构型并不是单学科气动特性最优的构型，而是气动特性和结构特性之间的一种折衷。

4.2 考虑静气动弹性影响的机翼剖面翼型优化设计

选取4.1节中的初始构型作为本节优化算例的初始构型，进行考虑静气动弹性影响的机翼剖面翼型优化设计研究，其中扭转角分布不变。如图 17所示，沿着机翼选取9个控制剖面，利用FFD方法，每个剖面有12个设计变量，结构设计仅仅进行承力结构的尺寸优化，设计变量描述如表 3所示。优化目标为航程最大，设计约束包括油箱容积不减，应力屈服极限(考虑2.5倍的过载系数)约束以及翼梢弯曲和扭转变形约束。为了对比研究考虑静气动弹性影响之后对机翼剖面翼型以及结构设计产生的影响，进行了以剖面翼型为设计变量的单学科的气动优化设计。最终，利用RANS方程对初始构型(Ori_aeroelasicity)、考虑静气动弹性影响的设计结果(Opt_aeroelasicity)以及单学科气动设计优化结果(Opt_aerodynamic)进行校核、对比分析。

如表 4所示，经过优化航程提高了3.05%，其中升阻比提高了1.53%，机翼重量减少了3.56%。优化后的构型翼梢挠度为1.28m(图 18)。从表面压力云图(图 19)和载荷分布(图 20)对比可以看出，考虑静气动弹性影响的设计结果相比于初始构型载荷适当向中段翼移动，内翼段和翼稍附近载荷相对较小。载荷分布的这种变化有利于在控制机翼翼根弯矩的情况下减小诱导阻力。相比之下，单学科的气动优化设计结果则表现为载荷的明显外移。

机翼剖面翼型的厚度、弯度、前加载和后加载对剖面当地载荷有着重要的影响。图 21显示考虑静气动弹性影响的设计结果相比于初始构型，中段翼最大厚度分布明显较小，同时图 23显示该处上翼面弯度基本没有减小，下翼面具有一定的前加载，而翼稍以及内翼段呈现出相反的趋势，因此载荷向中段翼移动。剖面压力分布对比图(图 24)显示，通过减小下翼面的厚度以及采用前加载或者后加载，可以使下翼面产生更多的升力，进而避免上翼面提供过大的升力，以此控制激波的强度。单学科的气动优化设计结果在外翼段翼型厚度减小，出现了明显的前加载和后加载, 上翼面激波强度被显著削弱。机翼翼盒面积分布(图 22)显示，为了满足油箱容积约束，无论是多学科的优化设计还是单学科的气动优化设计结果，翼根附近的翼盒面积都有所增大，以此弥补其它剖面翼盒面积的减小造成的油箱容积的损失，并且可以适当减小内翼段的载荷，使载荷外移。

从图 13和图 25可以看出，相比于初始构型，多学科优化结果除了kink附近蒙皮厚度有所增大外，整个外翼段蒙皮以及后梁厚度都明显减小。这主要是由于多学科优化设计结果的整个外翼段的翼型相对厚度明显增大(如图 23)，使得机翼翼盒承载能力增强，因此可以适当地减小外翼段蒙皮以及后梁厚度。在中段翼，虽然剖面翼型最大厚度减小，在一定程度上削弱了翼盒承载能力，但是后梁附近的机翼厚度并没有明显减小甚至有所增大，因此中段翼可以以较小的结构重量代价来维持机翼的承载能力。

5 结论

本文兼顾计算效率和计算复杂度之间的矛盾，建立了考虑静气动弹性影响的气动/结构一体化设计方法，并针对高马赫数跨声速客机，进行了考虑转角分布、剖面翼型和结构尺寸厚度影响的机翼多学科优化设计研究，得到的结论主要如下：

1) 通过进行CRM标模计算并与实验数据进行对比，说明基于全速式方程的数值方法可以作为初始设计阶段的数值评估手段，以显著地提高整个设计阶段的设计效率。同时，对F6模型进行的静气动弹性验证算例显示，建立的松耦合静气动弹性分析方法是可信的，可以用于考虑静气动弹性影响的优化设计研究中。

2) 在以扭转角为设计变量的算例中，相比于初始构型，多学科优化设计结果航程提高了5.63%。这主要得益于外翼段负扭转的减小，使得载荷向外翼段移动，升阻特性得到明显改善。相比于单学科的气动优化设计，考虑静气动弹性影响的多学科优化设计有效地控制了沿展向的扭转角分布，避免载荷分布过分外移。最终，以2.56%的结构重量损失代价获得了6.53%的升阻比的提高。

3) 在以剖面翼型为设计变量的算例中，相比于初始构型，多学科优化设计结果航程提高了3.05%。这主要得益于在保证油箱容积约束的条件下，通过调整机翼厚度分布以及前、后加载，在避免机翼激波强度明显增大的情况下，将载荷向中段翼移动。在此基础上，保证后梁所处位置附近的翼型厚度不减，并适当增大外翼段的翼盒面积，使得结构重量减少了3.56%，同时使升阻比略有提高(提高了1.53%)。