受加工精度、装配误差的影响，叶轮机械叶片实际外形和设计外形存在偏差，这种外形偏差的产生具有随机性，是叶轮机械典型不确定性影响因素之一。20世纪70年代研究人员即发现：外形偏差虽小，但是对气动性能的影响不容忽略[1-2]。随着气动设计技术的不断发展，通过优化设计提升气动性能越来越困难；相对而言，外形偏差引起的叶片性能变化较为明显，其气动不确定性影响研究已经引起了高度重视[3-10]。
基于几何外形设计参数和工况参数零偏差假设的气动设计较为成熟，相应的气动设计优化(ADO)称为确定性气动设计优化(DADO)。随着工业界对不确定性影响的日益重视，DADO面临着严峻挑战：考虑不确定性因素对气动性能的影响时，DADO外形在设计空间内是否仍然最优？答案是否定的。由于不确定性因素的客观性和随机性，最优气动外形的确定将依赖于设计外形在不确定性扰动空间内的统计气动性能，主要取决于：平均气动性能是否最优；气动性能参数对不确定性因素的敏感度是否最低。若DADO外形对不确定性因素比较敏感，外形偏差将导致气动性能在不确定性扰动空间内发生显著变化，甚至远离DADO解，意味着DADO外形在扰动空间内缺乏稳定性。在这种情况下，忽略不确定性影响的DADO很难满足实际工程需求，有必要开展鲁棒性气动设计优化研究。
鲁棒性是指系统在一定的系统参数或输入参数扰动下，维持系统输出稳定性的特性。鲁棒性优化在飞行器气动、结构强度设计等领域已经实现[11-13]，如：Paiva等[12]采用Kriging模型开展了考虑攻角及飞行高度不确定性影响的机翼多学科鲁棒性和可靠性优化设计，Ryan等[13]采用包括多目标优化在内的3种不同优化方法开展了考虑飞行参数不确定性影响的高超声速飞行器鲁棒性优化。飞行器鲁棒性优化往往忽略气动外形的不确定性影响，其主要原因是：气动外形偏差无法提前预估；相对于飞行器几何尺寸，外形偏差的影响可以忽略[12]。与飞行器不同的是：叶轮机械叶片几何尺寸较小，外形偏差对气动、结构等性能参数的不确定性影响不容忽略。整体上，叶轮机械鲁棒性气动设计优化研究起步也较晚[14-19]：南安普顿大学的Keane[14]开展了考虑叶片外形偏差影响的鲁棒性气动设计优化研究，比较了不同气动不确定性量化(UQ)方法对鲁棒性气动设计优化的影响；剑桥大学的Ghisu等[15]采用一维程序开展了压气机鲁棒性气动设计优化研究，采用混沌模型对级参数扰动的不确定性影响进行量化；Vinogradov等[17]开展了考虑边界扰动的高压涡轮传热和气动鲁棒性优化设计研究；Reis[18]、Wang[16]、Ma[19]等开展了考虑外形偏差的叶轮机械鲁棒性气动设计优化研究，采用多项式混沌模型进行气动UQ。
稳健性设计是鲁棒性设计的一种，稳健性气动设计优化(RADO)重点关注不确定性因素对目标气动函数的影响。如前文所述：RADO的核心技术之一是确定统计气动性能，即：确定目标气动函数的统计均值和方差。统计均值用于衡量平均气动性能，方差用于衡量气动性能的稳健性。因此，高效、高精度的气动UQ方法将直接影响RADO的优化精度和优化效率。目前在气动UQ中得到应用的方法有试验测量[1-2]、蒙特卡洛模拟(MCS)[3]、代理模型方法[8, 15, 18-19]、灵敏度分析方法[7, 9-10, 20-21]。整体上，对于外形偏差小尺度气动不确定性问题，基于灵敏度分析的气动UQ由于具有高效率、高精度的特点，目前已成为叶轮机械气动不确定性研究领域的热点。目前已公开发表的RADO研究[15-19]多采用基于代理模型的气动UQ方法来计算目标气动函数的统计均值和方差，导致优化效率较低。
本文将重点介绍叶轮机械RADO的基本原理及实现方法，并对某型跨声速涡轮叶栅进行气动优化设计，验证所发展的RADO方法的有效性和可靠性。首先介绍RADO需要解决的关键技术；之后介绍高效、高精度气动UQ方法及RADO目标函数梯度计算方法；最后开展跨声速涡轮叶栅DADO和考虑外形设计参数不确定性影响的RADO研究。通过优化气动外形降低叶栅总压损失系数及其方差，分析DADO与RADO对气动稳健性的影响机理。
1 跨声速叶栅数值模拟本文采用跨声速涡轮叶栅HS1A[22]作为研究对象，图 1为其几何外形及计算网格示意图，主要几何参数如表 1所示。该叶栅试验测量的工况参数为：来流角46°，进口为标准大气的压力和温度，出口等熵马赫数为1.06。本文只研究几何外形设计参数不确定性对气动优化设计的影响，暂不考虑边界流动参数的不确定性。进口给定总温、总压及来流角，出口给定反压。
|
| 图 1 HS1A几何示意图及计算网格 Fig. 1 Configuration and computation grid for HS1A |
| 参数 | 数值 |
| 弦长/mm | 40.0 |
| 轴向弦长/mm | 36.98 |
| 叶型前缘(后缘)角/(°) | 50.5(59) |
| 叶型安装角/(°) | 25.1 |
| 几何进口角/(°) | 46 |
| 栅距/mm | 29.14 |
流场数值模拟采用自开发程序求解定常流动控制方程，采用中心格式有限体积法、多重网格加速技术等。为了验证所采用数值模拟程序的有效性和可靠性，研究中首先求解Navier-Stokes方程和Spalart-Allmaras湍流模型方程，采用单块H-型网格，并对叶栅流向、周向网格进行收敛性分析。表 2给出了不同规模的网格及相应的总压损失系数，总压损失系数定义为
| $ \xi = \frac{{{p_{0,1}} - {p_{0,2}}}}{{{p_{0,2}} - {p_2}}} = \frac{{{p_{0,1}} - {p_2}}}{{{p_{0,2}} - {p_2}}} - 1 $ | （1） |
| 网格 | 流向 | 周向 | y+ | ζ |
| NS1 | 128 | 64 | 2.1 | 0.092 9 |
| NS2/NS5 | 160 | 64 | 2.1 | 0.094 8 |
| NS3 | 192 | 64 | 2.2 | 0.095 6 |
| NS4 | 160 | 48 | 5.2 | 0.093 3 |
| NS6 | 160 | 80 | 0.9 | 0.095 4 |
| Eu4 | 160 | 32 | 0.073 8 | |
| Eu5 | 160 | 48 | 0.076 2 | |
| Eu6 | 160 | 64 | 0.077 1 |
式中：p0, 1、p0, 2、p2分别表示进口总压、出口总压及静压。该叶栅的进口总压及出口静压均给定，通过计算出口总压即可确定总压损失系数。
由表 2可知，网格NS2/NS5的总压损失系数基本收敛。图 2给出了叶片表面等熵马赫数分布，由图可知，NS2计算结果与试验结果吻合较好。表 2及图 2结果表明，本研究所采用的自开发程序在模拟该叶栅流动时具有较高的可信度。
|
| 图 2 叶片表面等熵马赫数分布 Fig. 2 Isentropic Mach number distributions on blade surface |
后续研究中将采用“直接-伴随方法”计算二阶灵敏度(详细可见3.2节)，考虑到基于黏性伴随方法计算二阶灵敏度难度较大，目前还无法实现，研究中将求解定常Euler方程及相应的伴随方程计算二阶灵敏度。表 2给出了Euler方程求解的3套网格及总压损失系数，网格Eu5的总压损失系数基本收敛。图 2也给出了网格Eu5的等熵马赫数分布，与试验结果及NS2/NS5计算结果相比，Eu5所确定的激波位置有一定的偏差，其主要原因是Euler方程无法考虑流动黏性影响。
虽然忽略黏性影响的数值模拟及气动优化设计会给实际应用带来一定的局限性，但是基于定常Euler方程的流场计算作为一种流场数值模拟工具，在本文的RADO方法研究、RADO优化设计平台发展中仍具有重要作用。后续的流场数值模拟及气动优化设计，均采用网格Eu5。
2 稳健性气动优化设计 2.1 RADO基本原理确定性气动外形设计优化的前提条件是忽略不确定性因素对气动性能的影响，无需对目标气动函数在扰动空间内进行统计分析。DADO的数学模型可以定义为
| $ \begin{array}{*{20}{l}} {{\rm{Min }}}&{I(\mathit{\boldsymbol{w}},\mathit{\boldsymbol{b}})}\\ {{\rm{Sub}}}&{{c_i}(\mathit{\boldsymbol{w}},\mathit{\boldsymbol{b}}) \le 0,i = 1,2, \cdots ,{n_{\rm{c}}}} \end{array} $ | （2） |
式中：I表示目标气动函数；w、b分别表示流场、设计参数；ci为约束条件；nc为约束条件个数。式(2)为带约束最小化确定性优化问题。
与DADO不同的是：RADO关注的不再是某个特定解，而是不确定性参数扰动空间内的统计解。气动性能由统计均值表示，气动稳健性由方差表示。气动稳健性设计要求气动性能参数对不确定性因素的敏感度最低，也即气动性能的方差最小。RADO的数学模型可以定义为
| $ \begin{array}{*{20}{l}} {{\rm{Min }}}&{{\mu _I},{\sigma _I}}\\ {{\rm{Sub}}}&{{c_i}(\mathit{\boldsymbol{w}},\mathit{\boldsymbol{b}}) \le 0,i = 1,2, \cdots ,{n_{\rm{c}}}} \end{array} $ | （3） |
式中：μI和σI分别表示I的统计均值、方差。
对比式(2)、式(3)可知：DADO和RADO的基本流程一致，主要差别在于是否需要在不确定性参数扰动空间内对气动参数进行统计分析。图 3为DADO与RADO的气动函数统计特性示意图：假设考虑不确定性影响的气动函数变化均满足高斯分布，横坐标x表示一般化气动函数，纵坐标表示概率密度函数(PDF)，μ和σ分别表示气动函数x的统计均值和统计方差。以最小化优化问题为例，虽然DADO的气动函数统计均值比RADO更低，但是DADO的方差高于RADO，在这种情况下，RADO在保证优化叶片气动稳健性方面更具有优势。
|
| 图 3 DADO与RADO气动函数统计特性 Fig. 3 Statistics of aerodynamic function for DADO and RADO |
由上述对比分析可知：RADO是典型的多目标优化设计问题，其目标是：改善气动函数均值，提高优化叶片的平均气动性能；降低气动函数对不确定性因素的敏感度，提高优化叶片的气动稳健性。
2.2 叶栅气动优化设计总压损失系数的统计均值和方差是本文RADO的优化目标，研究中将采用加权方法构造多目标优化的目标函数：
| $ I = {\mu _\xi } + \lambda {\sigma _\xi } $ | （4） |
式中：λ为总压损失系数统计方差的权重，为无量纲参数。研究中，平均总压损失系数的权重为1，并保持不变；通过调整λ来对比分析方差权重对RADO优化结果的影响。
优化中暂不考虑气动约束，只约束叶栅厚度分布。本文将优化叶片中弧线外形并保持厚度分布不变。从前缘至98%弦长，在中弧线上分布N个型函数，中弧线变化量为
| $ {\rm{ \mathsf{ δ} }} y(x) = \sum\limits_{j = 1}^N {{V_j}} {b_j}(x) $ | （5） |
式中：Vj为型函数权重，即设计参数；bj(x)为Hicks-Henne型函数[23]，定义为
| $ \left\{ {\begin{array}{*{20}{l}} {{b_j}(x) = {\rm{si}}{{\rm{n}}^4}(\pi {x^{{p_j}}})}\\ {{p_j} = {\rm{ln}}(0.5)/{\rm{ln}}({x_{{\rm{c}},j}})\quad j = 1,2, \cdots ,N} \end{array}} \right. $ | （6） |
式中：xc, j为型函数控制点，沿中弧线轴向均匀分布6个控制点。
本文的RADO研究重点考虑几何外形设计参数的不确定性影响，也即Vj的不确定性影响。研究中假设所有设计参数彼此独立，其变化量满足均值为零的高斯分布：
| $ \begin{array}{l} f(\Delta V) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{\sqrt {2\pi } {\sigma _V}}}{\rm{exp}}\left( { - \frac{{{{(\Delta V)}^2}}}{{2\sigma _V^2}}} \right)}&{\frac{{\Delta V}}{{{\sigma _V}}} \in [ - E,E]}\\ 0&{{\rm{ otherwise }}} \end{array}} \right. \end{array} $ | （7） |
式中：ΔV表示设计参数变化量；E为高斯分布的截断边界；σV表示设计参数公差。在研究中，假设各设计参数的公差相同。
2.3 RADO方法及实现RADO目标气动函数的统计需要大量流场计算，气动UQ将占据优化设计的绝大部分计算时间，高效高精度的气动UQ方法是RADO需要重点关注的核心问题之一。在这种情况下，相对于优化设计方法的高效性，优化设计方法的精确性更值得关注。相对于气动参数的平均值，气动参数变化量较小，为了精确反映气动参数变化对设计优化的影响，必须采用高精度优化方法。目前在DADO中得到广泛应用的优化算法有：进化算法、代理模型、梯度方法。进化算法优化精度高，计算需求大；代理模型优化效率高，但是当模型固有的函数响应误差较大时，优化设计的信噪比较低，直接影响RADO的精度。梯度方法的精度及优化效率均较高。本文将采用原理简单、易实现的最速下降法进行优化设计。
作者的前期研究[9-10]表明，对于叶片几何外形偏差的不确定性问题，基于灵敏度分析的气动UQ具有高效、高精度的优势。本文的RADO研究仍将采用基于灵敏度分析的气动UQ方法来统计叶栅总压损失系数的均值和方差。基于一阶灵敏度分析的总压损失系数统计均值和方差分别为
| $ \left\{ {\begin{array}{*{20}{l}} {{\mu _\xi } = \frac{1}{m}\sum\limits_{i = 1}^m {({\xi _0} + \Delta {\xi _i})} = {\xi _0} + \frac{1}{m}\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^N {{g_j}} } \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \Delta {V_{ij}} = {\xi _0} + \sum\limits_{j = 1}^N {{g_j}} \left( {\frac{1}{m}\sum\limits_{i = 1}^m \Delta {V_{ij}}} \right) = {\xi _0}}\\ {\sigma _\xi ^2 = \frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {{\xi _0} + \sum\limits_{j = 1}^N {{g_j}} \Delta {V_{ij}} - {\mu _\xi }} \right)}^2}} = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{j = 1}^N {g_j^2} \left( {\frac{1}{m}\sum\limits_{i = 1}^m \Delta V_{ij}^2} \right) = \sigma _V^2\sum\limits_{j = 1}^N {g_j^2} } \end{array}} \right. $ | （8） |
式中：ξ0表示确定性总压损失系数；Δξ为总压损失系数变化；m为统计样本数；g表示总压损失系数对设计参数的一阶灵敏度。由式(8)可知：总压损失系数的统计均值与灵敏度无关，其方差与一阶灵敏度相关。计算确定所有设计参数的一阶灵敏度后，即可确定相应的方差。
通过气动UQ确定总压损失系数的统计均值和方差后，式(4)所示的目标函数即可计算确定。由式(8)可知，目标函数关于设计参数的梯度为
| $ {g_{I,i}} = \frac{{{\rm{ \mathsf{ δ} }} I}}{{{\rm{ \mathsf{ δ} }} {V_i}}} = \frac{{{\rm{ \mathsf{ δ} }} {\mu _\xi }}}{{{\rm{ \mathsf{ δ} }} {V_i}}} + \lambda \frac{{{\rm{ \mathsf{ δ} }} {\sigma _\xi }}}{{{\rm{ \mathsf{ δ} }} {V_i}}} = {g_i} + \lambda \frac{{\sigma _V^2}}{{{\sigma _\xi }}}\sum\limits_{j = 1}^N {{g_j}} {h_{ij}} $ | （9） |
式中：hij表示总压损失系数关于设计参数的二阶灵敏度；gI表示RADO目标函数关于设计参数的梯度。值得注意的是：在DADO中，由于σV=0，gI, i=gi。此外，由式(9)可知：当采用一阶灵敏度进行气动UQ时，RADO仍需要计算总压损失系数关于设计参数的二阶灵敏度。
当设计参数变化满足高斯分布且设计参数彼此独立(无关)时，基于二阶灵敏度分析的总压损失系数统计均值和方差为
| $ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {{\mu _\xi } = {\xi _0} + \frac{1}{m}\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^N {{g_j}} } \Delta {V_{ij}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{m}\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^N {\sum\limits_{k = 1}^N {\frac{1}{2}} } } {h_{jk}}\Delta {V_{ij}}\Delta {V_{ik}} = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\xi _0} + \frac{1}{2}\sum\limits_{j = 1}^N {{h_{jj}}} \left( {\frac{1}{m}\sum\limits_{i = 1}^m \Delta V_{ij}^2} \right) = {\xi _0} + \frac{1}{2}\sigma _V^2\sum\limits_{j = 1}^N {{h_{jj}}} } \end{array}\\ \sigma _\xi ^2 = \frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {\sum\limits_{j = 1}^N {{g_j}} \Delta {V_{ij}} - {\mu _\xi }} \right)}^2}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {\sum\limits_{j = 1}^N {\sum\limits_{k = 1}^N {\frac{1}{2}} } {h_{jk}}\Delta {V_{ij}}\Delta {V_{ik}} - {\mu _\xi }} \right)}^2}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma _V^2\sum\limits_{j = 1}^N {g_j^2} + \frac{1}{m}\sum\limits_{i = 1}^m {\left( {\sum\limits_{j = 1}^N {\sum\limits_{k = 1}^N {\frac{1}{2}} } {h_{jk}}\Delta {V_{ij}} \cdot } \right.} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\left. {\Delta {V_{ik}}} \right)^2} = \sigma _V^2\sum\limits_{j = 1}^N {g_j^2} + \frac{1}{2}\sigma _V^4\sum\limits_{j = 1}^N {\sum\limits_{k = 1}^N {h_{jk}^2} } \end{array} \right. $ | （10） |
由式(10)可知：总压损失系数的统计均值、方差和二阶灵敏度有关，计算式(4)所示的目标函数关于设计参数的梯度时，需要计算总压损失系数的三阶灵敏度。截至目前，二阶灵敏度的高效、高精度计算已经实现，三阶灵敏度计算还无法实现。因而，本文的RADO研究将采用基于一阶灵敏度的气动UQ方法，采用伴随方法计算总压损失系数的一阶、二阶灵敏度。
确定目标函数的梯度后，将采用最速下降法来实现优化设计。综上所述，基于梯度方法的叶栅RADO的流程图如图 4所示。图中：模块I为气动UQ，其核心是基于伴随方法的总压损失系数一阶、二阶灵敏度计算及基于一阶灵敏度分析的总压损失系数统计均值和方差计算；模块Ⅱ为优化设计模块。若忽略不确定性影响，模块I只保留流场、伴随场数值模拟及一阶伴随灵敏度计算等功能(无需进行基于二阶灵敏度的统计分析)，即可以实现基于梯度方法的叶栅DADO。
|
| 图 4 基于梯度方法的叶栅RADO流程 Fig. 4 Procedures of gradient-based RADO for cascades |
伴随方法由于灵敏度计算的高效、精确特点，在20世纪80年代由Jameson[24]应用于飞行器ADO研究。之后，该方法在飞行器ADO研究中得到了广泛应用[25-28]。21世纪初，该方法被应用于叶轮机械ADO研究[29]；近10年来，基于伴随方法的叶轮机械ADO得到了迅速发展[30-34]。接下来将简要介绍伴随方法的基本原理。
在计算气动函数关于设计参数的一阶灵敏度时，需要首先确定目标气动函数I、流动控制方程R关于设计参数b的一阶变分：
| $ {I = I(\mathit{\boldsymbol{w}},\mathit{\boldsymbol{b}}),\frac{{{\rm{ \mathsf{ δ} }} I}}{{{\rm{ \mathsf{ δ} }} {b_i}}} = \frac{{\partial I}}{{\partial \mathit{\boldsymbol{w}}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} + \frac{{\partial I}}{{\partial {b_i}}}} $ | （11） |
| $ {\mathit{\boldsymbol{R}}(\mathit{\boldsymbol{w}},\mathit{\boldsymbol{b}}) = {\bf{0}},\frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{R}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} = \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial \mathit{\boldsymbol{w}}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} + \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial {b_i}}} = {\bf{0}}} $ | （12） |
由式(11)可知：流场变分项δw/δbi是灵敏度计算的核心。
伴随方法基本原理是将流动控制方程作为约束引入到气动函数变分中，消除流场变分对灵敏度计算的影响。引入伴随算子Ψ，与式(12)相乘后再与式(11)相减，有
| $ \frac{{{\rm{ \mathsf{ δ} }} I}}{{{\rm{ \mathsf{ δ} }} {b_i}}} = \left( {\frac{{\partial I}}{{\partial \mathit{\boldsymbol{w}}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial \mathit{\boldsymbol{w}}}}} \right)\frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} + \left( {\frac{{\partial I}}{{\partial {b_i}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial {b_i}}}} \right) $ | （13） |
式中：若右端流场变分项的系数为零，则目标气动函数的一阶灵敏度与流场变分无关。同时需要求解关于伴随算子的线性方程组，此即为流动控制方程R的伴随方程：
| $ \frac{{\partial I}}{{\partial \mathit{\boldsymbol{w}}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial \mathit{\boldsymbol{w}}}} = 0 $ | （14） |
灵敏度计算为
| $ \frac{{{\rm{ \mathsf{ δ} }} I}}{{{\rm{ \mathsf{ δ} }} {b_i}}} = \frac{{\partial I}}{{\partial {b_i}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial {b_i}}} $ | （15） |
由式(15)可知：灵敏度计算只和设计参数变化引起的网格变分相关。相对于流场和伴随场计算，网格变分所需要的时间可以忽略。因而，基于伴随方法的一阶灵敏度计算所需要的时间约为两次流场计算时间，与设计参数个数基本无关。
3.2 二阶灵敏度计算直接-伴随方法[35-36]是目前最高效的二阶灵敏度计算方法，其基本原理介绍如下。
与一阶灵敏度计算类似，在进行二阶灵敏度计算时，需要首先确定目标气动函数I、流动控制方程R关于设计参数b的二阶变分。对式(11)、式(12)进行变分，有
| $ \begin{array}{l} \frac{{{{\rm{ \mathsf{ δ} }} ^2}I}}{{{\rm{ \mathsf{ δ} }} {b_i}{\rm{ \mathsf{ δ} }} {b_j}}} = \frac{{{\partial ^2}I}}{{\partial {b_i}\partial \mathit{\boldsymbol{w}}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_j}}} + \frac{{{\partial ^2}I}}{{\partial {b_j}\partial \mathit{\boldsymbol{w}}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} + \frac{{{\partial ^2}I}}{{\partial {b_i}\partial {b_j}}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{\partial ^2}I}}{{\partial {\mathit{\boldsymbol{w}}^2}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_j}}} + \frac{{\partial I}}{{\partial \mathit{\boldsymbol{w}}}} \cdot \frac{{{{\rm{ \mathsf{ δ} }} ^2}\mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}{\rm{ \mathsf{ δ} }} {b_j}}} \end{array} $ | （16） |
| $ \begin{array}{*{20}{c}} {\frac{{{{\rm{ \mathsf{ δ} }} ^2}\mathit{\boldsymbol{R}}}}{{{\rm{ \mathsf{ δ} }} {b_i}{\rm{ \mathsf{ δ} }} {b_j}}} = \frac{{{\partial ^2}\mathit{\boldsymbol{R}}}}{{\partial {b_i}\partial \mathit{\boldsymbol{w}}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_j}}} + \frac{{{\partial ^2}\mathit{\boldsymbol{R}}}}{{\partial {b_j}\partial w}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} + \frac{{{\partial ^2}\mathit{\boldsymbol{R}}}}{{\partial {b_i}\partial {b_j}}} + }\\ {\frac{{{\partial ^2}\mathit{\boldsymbol{R}}}}{{\partial {\mathit{\boldsymbol{w}}^2}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_j}}} + \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial \mathit{\boldsymbol{w}}}} \cdot \frac{{{{\rm{ \mathsf{ δ} }} ^2}\mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}{\rm{ \mathsf{ δ} }} {b_j}}} = \mathit{\boldsymbol{0}}} \end{array} $ | （17） |
根据伴随方法的一般原理，引入伴随算子Ψ，与式(17)相乘后再与式(16)相减，有
| $ \begin{array}{*{20}{l}} {{h_{ij}} = \mathit{\boldsymbol{A}}\frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} + \mathit{\boldsymbol{B}}\frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_j}}} + \mathit{\boldsymbol{C}}\frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}}} \cdot \frac{{{\rm{ \mathsf{ δ} }} \mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_j}}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{D}}\frac{{{{\rm{ \mathsf{ δ} }} ^2}\mathit{\boldsymbol{w}}}}{{{\rm{ \mathsf{ δ} }} {b_i}{\rm{ \mathsf{ δ} }} {b_j}}} + \mathit{\boldsymbol{E}}} \end{array} $ | （18） |
式中：hij表示二阶灵敏度；系数矩阵A、B、C、D及灵敏度矩阵E分别为
| $ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{A}} = \frac{{{\partial ^2}I}}{{\partial {b_j}\partial \mathit{\boldsymbol{w}}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{{\partial ^2}\mathit{\boldsymbol{R}}}}{{\partial {b_j}\partial \mathit{\boldsymbol{w}}}}}\\ {\mathit{\boldsymbol{B}} = \frac{{{\partial ^2}I}}{{\partial {b_i}\partial \mathit{\boldsymbol{w}}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{{\partial ^2}\mathit{\boldsymbol{R}}}}{{\partial {b_i}\partial \mathit{\boldsymbol{w}}}}}\\ {\mathit{\boldsymbol{C}} = \frac{{{\partial ^2}I}}{{\partial {\mathit{\boldsymbol{w}}^2}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{{\partial ^2}\mathit{\boldsymbol{R}}}}{{\partial {\mathit{\boldsymbol{w}}^2}}}}\\ {\mathit{\boldsymbol{D}} = \frac{{\partial I}}{{\partial \mathit{\boldsymbol{w}}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{\partial R}}{{\partial \mathit{\boldsymbol{w}}}}}\\ {\mathit{\boldsymbol{E}} = \frac{{{\partial ^2}I}}{{\partial {b_i}\partial {b_j}}} - {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\frac{{{\partial ^2}\mathit{\boldsymbol{R}}}}{{\partial {b_i}\partial {b_j}}}} \end{array}} \right. $ | （19） |
由式(18)可知：二阶灵敏度计算与流场的一阶、二阶变分有关；由式(19)可知：系数矩阵D为伴随方程，即D=0。因而，流场二阶变分对二阶灵敏度计算没有影响。
由上述分析可知：只需要N+1次流场计算确定对应每个设计参数的流场变分，N为设计参数个数；再进行一次伴随方程计算，即可以根据式(18)计算所有二阶灵敏度。伴随方程与流动控制方程数值求解所需要的计算时间基本相当，因而，采用上述直接-伴随方法，二阶灵敏度计算所需要的总流场计算次数约为N+2。此外，由于N+1次流场计算结果已经确定，可以采用直接差分法(FDM)计算一阶灵敏度。
基于直接-伴随方法的二阶灵敏度计算精度在前期工作中已经与基于FDM的二阶灵敏度进行了对比，其精度是可靠的[10]。此外，基于二阶灵敏度分析的气动UQ结果与MCS结果进行了对比，气动参数变化量的统计结果精度也得到了验证。因此，直接-伴随方法将用于式(9)所示的RADO目标函数梯度的计算。基于Euler方程的直接-伴随方法详细推导及二阶灵敏度计算请参考本文作者前期工作[10]。
4 优化设计为了对比分析RADO在改善优化叶片气动稳健性方面的特点，首先开展基于梯度方法的HS1A叶栅DADO研究。大量不确定性研究表明[3, 7-10]：几何外形变化公差对叶片气动性能不确定性变化有重要影响。因而，研究中将考虑不同的设计参数变化公差对优化设计的影响。此外，由前文分析可知：RADO是典型的多目标优化设计问题，不同优化目标间的权重对优化设计结果也有重要影响。如式(4)所示，研究中将考虑不同的方差权重λ对RADO的影响。表 3给出了4种不同类型的RADO。
| 参数 | C1 | C2 | C3 | C4 |
| σV/mm | 0.025 | 0.050 | 0.025 | 0.050 |
| λ | 1.0 | 1.0 | 5.0 | 5.0 |
DADO与RADO的本质区别在于是否考虑几何设计参数不确定性对总压损失系数变化的影响。这种影响将反映在优化设计目标函数对设计参数的梯度上。图 5对比了DADO与RADO第1步优化设计目标函数的梯度，图中横坐标V表示设计参数，纵坐标G表示梯度。由图可知：前5个设计参数的梯度几乎完全重合，第6个设计参数(最靠近尾缘)的梯度有非常明显的差异。
|
| 图 5 第1步优化设计梯度对比 Fig. 5 Comparisons of gradients in first optimization cycle |
由式(9)可知：导致目标函数梯度差异的主要原因有：设计参数公差、方差权重、总压损失系数的一阶、二阶灵敏度。对于不同的设计参数，当公差和气动函数的方差系数确定时，导致RADO与DADO梯度差异的主要因素是总压损失系数对设计参数的敏感度。第6个设计参数处于跨声速流动区域，气动外形的微小变化即可带来总压损失系数的较大变化。此外，对比不同的RADO梯度可知：随着方差权重的增加，目标函数梯度增加；随着几何设计参数公差的增加，目标函数梯度也增加。整体上，RADO梯度与DADO梯度的差异大小，反映了设计参数不确定性变化对气动函数变化的影响大小。
图 6是DADO与RADO的目标函数、总压损失系数收敛曲线。80步后优化基本收敛，收敛的目标函数几乎重合，而总压损失系数有明显差异：RADO优化总压损失系数比DADO低，且随着公差和方差权重的增大，总压损失系数下降更多。当λ=5时，RADO总压损失系数明显低于DADO。
|
| 图 6 优化设计收敛曲线 Fig. 6 Convergence history of optimization design |
由RADO的基本原理可知：叶片气动外形优化设计除了改善平均气动性能，还要降低气动参数对设计参数变化的敏感度，也即降低总压损失系数的方差。采用基于一阶灵敏度分析的气动UQ方法，总压损失系数的统计方差与一阶灵敏度相关，如式(8)所示。由于总压损失系数对设计参数的一阶灵敏度已知，DADO中也可以计算总压损失系数的方差。图 7给出了DADO与RADO的总压损失系数统计方差收敛曲线，其中：DADO_C1和DADO_C2分别表示公差为0.025 mm和0.050 mm的统计结果。由图可知：随着设计参数变化的公差减小，总压损失系数的方差降低；当设计参数公差相同时，随着方差权重的增加，总压损失系数方差下降；RADO的优化方差明显低于DADO优化方差，反映了RADO在改善优化叶片气动稳健性方面的优势。
|
| 图 7 总压损失系数方差收敛曲线 Fig. 7 Convergence history of standard deviation of total pressure loss coefficients |
表 4给出了DADO与RADO的主要优化结果，包括：平均总压损失系数、总压损失系数统计方差。与相同公差的DADO相比，RADO的总压损失系数统计均值和方差均明显下降，尤其是方差下降更明显。随着方差权重系数λ的增加，统计均值下降更多，意味着总压损失系数偏差更大；同时，统计方差下降也更多，意味着总压损失系数对设计参数变化的敏感度降低，证实了稳健性优化设计的有效性。
| 优化方法 | μ/10-2 | σ/10-4 | Δμ/% | Δσ/% |
| DADO_C1 | 5.748 | 1.564 | ||
| DADO_C2 | 5.748 | 3.128 | ||
| RADO_C1 | 5.737 | 1.498 | -0.19 | -4.22 |
| RADO_C2 | 5.729 | 2.910 | -0.33 | -6.97 |
| RADO_C3 | 5.703 | 1.270 | -0.78 | -18.8 |
| RADO_C4 | 5.671 | 2.077 | -1.34 | -33.6 |
图 8给出了优化叶片与原始叶片的气动外形及叶面等熵马赫数分布。优化后中弧线在叶片中部的曲率明显减小，抑制流动加速，如图 8(b)所示：叶片中部吸力面等熵马赫数明显降低，从而降低了激波前马赫数；同时，位于原始叶片吸力面65%弦长处的第1道激波完全消失。此外，优化后叶片后半部分中弧线曲率明显增加，流动进一步加速，如图 8(b)所示，优化叶片80%弦长后仍是超声速区，位于原始叶片吸力面80%弦长处的第2道激波强度显著减弱，有利于降低总压损失系数。整体上，随着设计参数公差和方差权重的增加，上述中弧线曲率、吸力面马赫数的变化趋势更加明显。
|
| 图 8 优化对比 Fig. 8 Optimization comparisons |
由上述结果可知：RADO优化叶片的总压损失系数下降，方差也低于DADO，气动稳健性得到了明显改善。引起方差下降的主要原因是：外形优化设计降低了激波强度，降低了总压损失系数对外形变化的敏感度。相对于DADO，在相同的公差条件下，RADO总压损失系数方差下降更多，主要原因在于激波强度的减弱。
为了进一步揭示RADO在改善叶片气动稳健性方面的优势，接下来对原始叶片、DADO、RADO_C2和RADO_C4优化叶片的流场进行统计分析，几何设计参数的公差均为0.05 mm。叶片样本数为1 000，采用相同的流场数值模拟程序计算叶片总压损失系数，之后再进行MCS统计分析。定义总压损失系数相对变化量为
| $ {\rm{ \mathsf{ δ} }} \xi = \left( {1 - \frac{\xi }{{{\xi _0}}}} \right) \times 100\% $ | （20） |
式中：ξ0表示4种基准叶片的总压损失系数。
图 9对比了不同叶片的总压损失系数相对变化量的统计值，图中横坐标1~4分别对应上述4种叶片。由图可知：相对于原始叶片，气动外形优化设计使总压损失系数变化量的统计值显著下降；此外，相对于DADO，RADO的统计均值略高，RADO的统计方差较低；随着方差权重的增加，RADO统计均值升高，统计方差下降，这是由多目标优化设计中各目标的竞争关系导致的。整体上，通过MCS统计分析进一步证实了RADO在改善气动稳健性方面的优势。
|
| 图 9 不同叶片的总压损失系数变化量统计值 Fig. 9 Statistics of total pressure loss coefficient change for different blades |
为了深入分析RADO对叶片气动稳健性的影响机理，接下来还将对比上述4种叶片的统计流场。图 10、图 11分别给出了上述4种叶片的等熵马赫数变化的统计均值、方差云图。与原始叶片相比，DADO与RADO均显著降低了吸力面激波附近的等熵马赫数变化统计值。统计均值与统计方差的降低均表明：叶片气动外形的变化对激波附近的流动影响下降，流场对设计参数变化的敏感度降低。然而，由图 10、图 11的云图很难发现RADO在改善等熵马赫数变化量统计值方面的优势。
|
| 图 10 通道等熵马赫数变化均值云图 Fig. 10 Contours of mean change of isentropic Mach number in blade passage |
|
| 图 11 通道等熵马赫数变化方差云图 Fig. 11 Contours of standard deviation of isentropic Mach number in blade passage |
图 12、图 13分别为DADO与RADO_C2 (σV=0.05 mm, λ=1.0)、RADO_C4(σV=0.05 mm, λ=5.0)优化叶片的等熵马赫数变化量的统计均值和方差(下文统称为统计量)在叶片表面的分布。图中：SS、PS分别表示叶片吸力面、压力面。为了更好地对比分析DADO与RADO在提高气动稳健性方面的特点，图 12(b)和图 13(b)给出了引入放大因子后的叶片吸力面统计量分布。放大因子主要用于调整RADO与DADO统计量之间的偏差，即
| $ {f_{{\rm{RADO}}}} = {f_{{\rm{DADO}}}} + a({f_{{\rm{0,RADO}}}} - {f_{{\rm{DADO}}}}) $ | （21） |
|
| 图 12 叶面等熵马赫数变化统计均值分布 Fig. 12 Distributions of means of isentropic Mach number change on blade surface |
|
| 图 13 叶面等熵马赫数变化统计方差分布 Fig. 13 Distributions of variances of isentropic Mach number change on blade surface |
式中：f表示统计量；a为放大因子；下标0表示原始统计量。
由图 12(a)和图 13(a)可知：相对于原始叶片，优化叶片的等熵马赫数变化量的统计均值和方差显著降低，表明：DADO与RADO均能降低流场对设计参数变化的敏感性；在压力面上，DADO与RADO的统计量非常接近，主要是因为压力面流动受激波影响较小，其流场对设计参数变化欠敏感。图 12(b)和图 13(b)中，左图对比了激波区域(60%~85%轴向弦长)的统计分布量，右图对比了尾缘附近(85%~100%轴向弦长)的统计量分布。由图 12(b)和图 13(b)可知：在激波区域RADO的统计量明显低于DADO，且随着方差权重的增加，RADO统计量下降更明显。可以预见的是：随着方差权重的增加，激波附近等熵马赫数变化量的统计均值和方差将继续下降，流场对设计参数变化的敏感性将进一步降低。在尾缘附近，RADO的均值比DADO高；RADO的方差在峰值前明显低于DADO，在峰值后略高于DADO。导致这种结果的原因可能是涡轮叶栅厚尾缘附近流场数值模拟精度欠佳。
由上述统计分析可知：DADO及RADO由于减弱了激波强度，使得流场对设计参数变化的敏感度下降，从而使得优化叶片的气动稳健性得到改善。在无黏流动中，流动损失的主要来源是激波。由于DADO也能降低激波强度，激波附近流场对设计参数变化的敏感度降低，导致总压损失系数的统计方差也随之降低，改善了优化叶片的气动稳健性。此外，由式(8)可知：总压损失系数的统计方差与一阶灵敏度及设计参数公差相关；在给定的设计参数公差条件下，随着优化设计的进行，总压损失系数的一阶灵敏度不断下降，也将导致其统计方差的降低及气动稳健性的改善。相对于DADO，RADO由于显式地考虑了总压损失系数统计方差的优化，能够更有效地降低总压损失系数的统计方差，进一步提高优化叶片的气动稳健性。
由于目前的研究只考虑了激波强度变化对叶片气动稳健性的影响，DADO与RADO改善平均气动性能及气动稳健性的趋势是一致的。在更复杂的流动中，如同时考虑激波、边界层及激波边界层干涉等，DADO与RADO对气动稳健性的影响机理有待进一步深入分析。
5 结论本文介绍了叶轮机械RADO基本原理，并针对RADO的核心问题：气动UQ及优化方法进行分析。采用伴随方法计算目标气动函数关于设计参数的一阶、二阶灵敏度，高效、高精度地实现了目标函数关于设计参数的梯度计算，实现了考虑几何设计参数不确定性影响的稳健性优化设计。开展了某型跨声速叶栅的DADO及RADO研究，并进行了详细对比分析，主要结论有：
1) 在无黏流动中，DADO与RADO均能降低跨声速叶栅的总压损失系数，并降低总压损失系数的方差，改善优化叶片的气动稳健性；相对于DADO，RADO能够更有效地降低总压损失系数的统计方差，其优化叶片气动稳健性更优。
2) 随着设计参数公差的增加，总压损失系数方差也增加；随着统计均值、方差多目标优化设计中方差权重的增加，总压损失系数下降更明显，优化叶片的气动稳健性更优。
3) 无黏流动中，影响叶片气动稳健性的主要因素是激波；气动外形优化设计通过减弱激波强度降低流场对设计参数变化的敏感度，从而降低跨声速叶栅总压损失系数变化量的统计值，改善优化叶片的稳健性。
| [1] | BAMMERT K, STOBBE H. Results of experiments for determining the influence of blade profile changes and manufacturing tolerances on the efficiency, the enthalpy drop, and the mass flow of multi-stage axial turbines: ASME Paper No. 70 WA/GT-4[R]. New York: ASME, 1970. |
| [2] | BAMMERT K, SANDSTEDE H. Influence of manufacturing tolerances and surface roughness of blades on the performance of turbines[J]. Journal of Engineering for Power, 1976, 98(1): 29-36. |
| Click to display the text | |
| [3] | GARZON V, DARMOFAL D. Impact of geometric variability on axial compressor performance[J]. Journal of Turbomachinery, 2003, 125(4): 692-703. |
| Click to display the text | |
| [4] | EDWARDS R, ASGHAR A, WOODASON R, et al. Numerical investigation of the influence of real world blade profile variation on the aerodynamic performance of transonic nozzle guide vanes[J]. Journal of Turbomachinery, 2012, 134(2): 021014. |
| Click to display the text | |
| [5] | LANGE A, VOIGT M, VOGELER K, et al. Impact of manufacturing variability on multistage high-pressure compressor performance[J]. Journal of Engineering for Gas Turbines and Power, 2012, 134(11): 112601. |
| Click to display the text | |
| [6] | SCHNELL R, LENGYEL-KAMPMANN T, NICKE E. On the impact of geometric variability on fan aerodynamic performance, unsteady blade row interaction, and its mechanical characteristics[J]. Journal of Turbomachinery, 2014, 136(9): 091005. |
| Click to display the text | |
| [7] | YANG J, XIONG J, MCBEAN I, et al. Performance impact of manufacturing variations for multi-stage steam turbines[J]. Journal of Propulsion and Power, 2017, 33(4): 1031-1036. |
| Click to display the text | |
| [8] |
蔡宇桐, 高丽敏, 马驰, 等. 基于NIPC的压气机叶片加工误差不确定性分析[J]. 工程热物理学报, 2017, 38(3): 490-497. CAI Y T, GAO L M, MA C, et al. Uncertainty quantification on compressor blade considering manufacturing error based on NIPC method[J]. Journal of Engineering Thermophysics, 2017, 38(3): 490-497. (in Chinese) |
| Cited By in Cnki (9) | Click to display the text | |
| [9] |
罗佳奇, 朱亚路, 刘锋. 基于伴随方法的叶片加工偏差气动灵敏度分析[J]. 工程热物理学报, 2017, 38(3): 498-504. LUO J Q, ZHU Y L, LIU F. Aerodynamic sensitivity analysis for manufacturing variations of a turbine blade by an adjoint method[J]. Journal of Engineering Thermophysics, 2017, 38(3): 498-504. (in Chinese) |
| Cited By in Cnki (2) | Click to display the text | |
| [10] | LUO J, LIU F. Statistical evaluation of performance impact of manufacturing variability by an adjoint method[J]. Aerospace Science and Technology, 2018, 77: 471-484. |
| Click to display the text | |
| [11] | LI H, MA C. Hybrid dimension-reduction method for robust design optimization[J]. AIAA Journal, 2013, 51(1): 138-144. |
| Click to display the text | |
| [12] | PAIVA R, CRAWFORD C, SULEMAN A. Robust and reliability-based design optimization framework for wing design[J]. AIAA Journal, 2014, 52(4): 711-724. |
| Click to display the text | |
| [13] | RYAN K, LEWIS M, YU K. Comparison of robust optimization methods applied to hypersonic vehicle design[J]. Journal of Aircraft, 2015, 52(5): 1510-1523. |
| Click to display the text | |
| [14] | KEANE A J. Comparison of several optimization strategies for robust turbine blade design[J]. Journal of Propulsion and Power, 2009, 25(5): 1092-1099. |
| Click to display the text | |
| [15] | GHISU T, PARKS G, JARRETT J, et al. Robust design optimization of gas turbine compression systems[J]. Journal of Propulsion and Power, 2011, 27(2): 282-295. |
| Click to display the text | |
| [16] | WANG X, HIRSCH C, LIU Z, et al. Uncertainty-based robust aerodynamic optimization of rotor blades[J]. International Journal for Numerical Methods in Engineering, 2013, 94: 111-127. |
| Click to display the text | |
| [17] | VINOGRADOV K, KRETININ G, OTRYAHINA K, et al. Robust optimization of the hpt blade cooling and aerodynamic efficiency: ASME Paper No. GT2016-56195[R]. New York: ASME, 2016. |
| [18] | REIS C, MANZANARES-FILHO N, DE LIMA A. Robust optimization of turbomachinery cascades using inverse methods[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2016, 38(1): 297-305. |
| Click to display the text | |
| [19] | MA C, GAO L, CAI Y, et al. Robust optimization design of compressor blade considering machining error: ASME Paper GT2017-63157[R]. New York: ASME, 2017. |
| [20] | GIEBMANNS A, BACKHAUS J, FREY C, et al. Compressor leading edge sensitivities and analysis with an adjoint flow solver: ASME Paper GT2013-94427[R]. New York: ASME, 2013. |
| [21] | ZAMBONI G, BANKS G, BATHER S. Gradient-based adjoint and design of experiment CFD methodologies to improve the manufacturability of high pressure turbine blades: ASME Paper GT2016-56042[R]. New York: ASME, 2016. |
| [22] | JOUINI D. Experimental investigation of two transonic linear turbine cascades at off-design conditions[D]. Ottawa: Carleton University, 2000. |
| [23] | HICKS R, HENNE P. Wing design by numerical optimization[J]. Journal of Aircraft, 1978, 15(7): 407-412. |
| Click to display the text | |
| [24] | JAMESON A. Aerodynamic design via control theory[J]. Journal of Scientific Computing, 1988, 3(3): 233-260. |
| Click to display the text | |
| [25] | NADARAJAH S, JAMESON A. Optimum shape design for unsteady flows with time-accurate continuous and discrete adjoint methods[J]. AIAA Journal, 2007, 45(7): 1478-1491. |
| Click to display the text | |
| [26] |
熊俊涛, 乔志德, 杨旭东, 等. 基于黏性伴随方法的跨声速机翼气动优化设计[J]. 航空学报, 2007, 28(2): 281-285. XIONG J T, QIAO Z D, YANG X D, et al. Optimum aerodynamic design of transonic wing based on viscous adjoint method[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(2): 281-285. (in Chinese) |
| Cited By in Cnki (31) | Click to display the text | |
| [27] | MARTINS J, LAMBE A. Multidisciplinary design optimization:A survey of architectures[J]. AIAA Journal, 2013, 51(9): 2049-2075. |
| Click to display the text | |
| [28] |
黄江涛, 周铸, 刘刚, 等. 飞行器气动/结构多学科延迟耦合伴随系统数值研究[J]. 航空学报, 2018, 39(5): 121731. HUANG J T, ZHOU Z, LIU G, et al. Numerical study of aero-structural multidisciplinary lagged coupled adjoint system for aircraft[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(5): 121731. (in Chinese) |
| Cited By in Cnki (13) | Click to display the text | |
| [29] | YANG S, WU H, LIU F, et al. Aerodynamic design of cascades by using an adjoint equation method: AIAA-2003-1068[R]. Reston: AIAA, 2003. |
| [30] | WANG D, HE L. Adjoint aerodynamic design optimization for blades in multi-stage turbo-machines:part i-methodology and verification[J]. Journal of Turbomachinery, 2010, 132(2): 021011. |
| Click to display the text | |
| [31] | LUO J, XIONG J, LIU F, et al. Three-dimensional aerodynamic design optimization of a turbine blade by using an adjoint method[J]. Journal of Turbomachinery, 2011, 133(1): 011026. |
| Click to display the text | |
| [32] | WALTHER B, NADARAJAH S. Constrained adjoint-based aerodynamic shape optimization of a single-stage transonic compressor[J]. Journal of Turbomachinery, 2013, 135(2): 021017. |
| Click to display the text | |
| [33] | LUO J, LIU F, MCBEAN I. Turbine blade row optimization through endwall contouring by an adjoint method[J]. Journal of Propulsion and Power, 2015, 31(2): 505-518. |
| Click to display the text | |
| [34] | MA C, SU X, YUAN X. An efficient unsteady adjoint optimization system for multistage turbomachinery[J]. Journal of Turbomachinery, 2016, 139(1): 011003. |
| Click to display the text | |
| [35] | PAPADIMITRIOU D, GIANNAKOGLOU K. Comput-ation of Hessian matrix in aerodynamic inverse design using continuous adjoint formulations[J]. Computers & Fluids, 2008, 37: 1029-1039. |
| Click to display the text | |
| [36] | PAPADIMITRIOU D, GIANNAKOGLOU K. The continuous adjoint approach for second order sensitivities in viscous aerodynamic inverse design problems[J]. Computers & Fluids, 2009, 38: 1539-1548. |
| Click to display the text |