基于数论网格法与Morris法的翼型气动特性敏感性分析

文泽军，孟祥恒，肖 钊，张 帆

（湖南科技大学机电工程学院，湖南 湘潭 411201）

翼型气动特性是风力机设计的决定性因素［1］。通过改变翼型的厚度、弯度和前缘半径等设计变量，可以设计出不同几何外形和气动性能的风力机翼型［2］。目前，CST（class function/shape function transformation，类型函数/形状函数转换）方法被广泛应用于风力机翼型设计中，提高了翼型设计的效率［3］。在CST方法的应用过程中，参数化模型的设计变量通常在一定区间内变化，且变化带有一定的随机性，常被认为是服从正态分布的随机变量［4］。该随机变量会导致所设计的风力机翼型出现几何偏差，从而可能导致风力机翼型的气动性能不能满足设计要求［5］。如Campobasso等［6］指出风力机翼型在设计及制造时随机产生的几何误差会对风力机的功率特性造成一定影响。因此，在风力机翼型气动设计中须考虑几何不确定性的影响。国内外学者围绕翼型气动特性开展了大量研究。如：Ju等［7］基于蒙特卡洛模拟（Monte Carlo simulation，MCS）研究了风力机翼型气动特性对翼型前缘的敏感性，但MCS法的样本收敛性易受到模拟次数及样本数的限制；
Ernst等［8］采用拉丁超立方法研究了风力机翼型的相对厚度和相对弯度等几何结构参数对气动特性的影响，但拉丁超立方法存在样本空间收敛性不足、不适用于高维运算等问题；
王晓东等［9-10］研究了多项式混沌方法在随机方腔流动模拟中的应用；
姜昌伟等［11］采用嵌入式多项式混沌展开法研究了在随机边界下流动与传热问题的不确定性，但由于须对系统内部方程进行扩维处理，容易引入数值误差；
刘全等［12］研究了非嵌入式多项式混沌方法在拉氏计算中的应用，但在正交多项式基的构建中易出现矩匹配方程的奇异解，难以保证计算精度；
为了较好地求解圆筒实验模型，梁霄等［13］提出了一种多元非嵌入多项式混沌方法，但在预测非嵌入样本点处的响应值时，样本点的数量和位置都不是任意的，求解展开式系数时主观性较大，容易引入近似误差；
邬晓敬等［14］采用非嵌入多项式混沌方法开展翼型气动特性全局敏感性分析，但因易受多项式模型的阶数影响，近似解计算的主观性较大；
祝宏宇等［15］采用Sobol全局敏感性分析方法研究翼型气动特性的敏感性，但该方法在求解高维运算问题时时间成本过高；
Wang等［16］采用多点泰勒展开的近似伴随方法分析翼型几何的设计变量对气动特性的影响，但近似伴随方法依然易受到多点泰勒展开式项数的影响，导致样本空间的计算精度不高；
Peter等［17］对翼型气动设计中的敏感性分析方法进行了总结。综上可知，目前的研究方法存在计算效率低、样本收敛性差及不适用于高维运算等不足，因此须对翼型气动特性的敏感性作进一步研究。

近年来，数论网格（number theoretic net，NT-net）法被广泛应用于工程统计分析领域。方开泰等［18］阐述了NT-net法在统计中的应用；
刘飞等［19］阐述了NT-net法在极大似然估计中的应用，指出NT-net法具有较高的采样效率和计算精度；
赖斌等［20］阐述了NT-net法在流体随机模拟中的应用，指出NT-net法可适用于高维随机模拟。目前，将NT-net法应用于风力机翼型气动特性敏感性分析的研究还不多见。

因此，本文以风力机翼型S832（以下简称翼型）为研究对象，基于CST方法构建翼型参数化模型，并提出一种基于NT-net法与Morris法的翼型气动特性敏感性分析方法。首先，基于CST方法建立精度较高的翼型参数化模型；
其次，应用NT-net法对参数化模型中的设计变量进行抽样，并采用改进的Morris法分析NT-net法的抽样结果，研究在翼型设计过程中各设计变量对翼型气动特性的贡献程度，并通过MCS进行验证分析；
最后，研究主要设计变量对翼型外形及气动特性的影响，以期为翼型的设计提供理论参考。

1.1 翼型参数化模型

1.1.1 翼型的几何特征

翼型为风力机叶片沿展向的某一横截面的形状，其前缘呈圆弧形，后缘呈钝形或尖形。弦线与前缘的交点为前缘点，弦线与后缘的交点为后缘点。其几何特征与气动特性密切相关。翼型的几何特征如图1所示。

图1 翼型几何特征示意Fig.1 Schematic of geometric characteristics of airfoil

1.1.2 翼型参数化方法

CST方法是根据翼型的几何特征，通过类函数和型函数来表示翼型外形变化的参数化方法。其优势是可以用少量的变量表达出可由计算机识别的翼型几何形状。因此，本文采用CST方法建立翼型参数化模型。以前缘点为原点O，以前缘点指向后缘点的方向为横向即ψ向的正向，以垂直于弦线向上为纵向即ξ向的正向，建立坐标系O-ψξ，则翼面可以表示为：

式中：B(ψ )为类函数；
G(ψ )为型函数；
ΔLe为翼型后缘厚度；
c为翼型的弦长。

其中：

式中：N1、N2为类函数系数，常取N1=0.5，N2=1.0；
t为翼型Bernstein多项式阶数，t=0，1，…，W；
At+1为Bernstein多项式的待定系数，可由最小二乘法拟合得到。

参数化误差Δy为：

式中：yψ为参数化翼型外形中横坐标为ψ处对应的纵坐标。

1.2 翼型气动特性

翼型气动特性主要由升力系数Cl和阻力系数Cd表征。升阻比ε定义为升力系数与阻力系数的比，即：

Xfoil软件是一款交互式独立翼型的设计与分析软件。与其他CFD（computational fluid dynamics，计算流体动力学）模拟软件相比，Xfoil软件具有收敛快、效率高和易操作等优点。基于线性变换的高阶面元法，可利用Xfoil软件计算Cl和Cd等。

2.1 NT-net法

数论方法（number theory method，NTM）的实质是在b维单位立方体Mb上找到一个均匀散布的点集［21］，该点集通常称为NT-net，偏差较小的点集称为好格子点（good lattice point，GLP）集合。

设 P={gk(k=1，…，n)}为 Mb上的一个点集。对于r∈Mb，令U(r，P)为P中满足gk≤r的点数，则P的偏差D（n，P）为：

式中：|0，r|为边长为r的矩形，R(|0，r|)为矩形|0，r|的体积空间。

确定最优系数(h1，…，hb)是获得GLP集合的关键环节。当维数b≤18时，应用生成矢量产生GLP集合。

令生成矢量 (n；
h1，…，hb)满足最大公约数(n，hi)=1，i=1，…，b，1≤hi

令 qki满 足 1≤qki≤n，则 具 有 最 小 偏 差 的 集 合{gk=(gk1，… ，gkb)} 即为(n；
h1，···，hb)的GLP集合。

2.2 基于NT-net法抽样的敏感性分析方法

抽样是指将b维立方体Mb的GLP集合的点映射到给定分布上的过程。敏感性分析常通过MCS来抽样模拟各种输入参数的可变性［22］。NT-net法被称为伪MCS法，因其良好的样本收敛性而被广泛应用。Morris法是由Morris提出的全局敏感性分析方法，其敏感性指标为变量元效应的均值和标准差，可用来评价系统输入量对系统输出量的影响［23］。

为量化每个Bernstein多项式待定系数的不确定性影响，令待定系数为服从正态分布N(μ，σ)的设计变量 Z，Z=(z1，z2，…，zb)T，Z 的联合分布函数为F(zi)。根据样本数n及维数b，通过查找生成矢量或采用分圆域方法得到b维立方体Mb的NT-net及GLP集合，则GLP集合中的每个点gki表征b维设计变量Z的一个样本值，由此可以建立翼型参数化模型的设计空间。

建立敏感性模型 Q=f(Z)，其中 zi～N(μ，σ )。对于通过NT-net法抽样得到的gki，在第v次抽样的变量元效应Ei，v及其均值 μi与标准差Si为：

针对元效应会在非单调函数中同时出现正负值的问题，Campolongo等［24］将均值μi改进为新均值：

敏感性分析的首要目标是计算并选择均值较高或标准差较低的变量元效应。因此，定义变异系数ζ为新的敏感性指标，其为变量元效应的均值与标准差的比值，即：

2.3 基于NT-net法抽样的敏感性分析流程

翼型气动特性敏感性分析的流程如图2所示，具体为：

图2 翼型气动特性敏感性分析流程Fig.2 Sensitivity analysis flow of aerodynamic characteristics of airfoil

1）确定设计变量Z、样本数n和维数b，确定设计变量的分布特征；

2）根据维数b的大小选定适用方法进行NT-net抽样并生成对应设计变量Z的GLP集合；

3）从设计变量Z的GLP集合中抽取样本点，根据式（1）至式（5）构建翼型参数化模型的设计空间，利用Xfoil软件求解翼型参数化模型的气动特性；

4）基于改进的Morris法，分析设计变量的变化对气动特性的影响，根据式（12）至（16）计算翼型气动特性的变异系数ζ，并找出影响翼型气动特性的主要设计变量。

3.1 翼型参数化模型的构建及气动特性敏感性分析

3.1.1 翼型参数化模型的构建

选取风力机专用翼型S832为研究对象。为了得到拟合精度较高的翼型参数化模型，取Bernstein多项式阶数W=4，5，6，采用CST方法分别对翼型的上、下翼面进行拟合并构建参数化模型。将翼型的几何外形坐标导入翼型设计软件Profili，可以更直观地表示翼型的几何外形特征。在Profili中，纵坐标y/c表示翼型几何外形的纵坐标y相对于翼型弦长c的位置关系，横坐标x/c表示翼型几何外形的横坐标x相对于翼型弦长c的位置关系。翼型模型的参数化误差如图3所示。

图3 翼型模型的参数化误差Fig.3 Parameterized error of airfoil model

由图3可知，Bernstein多项式阶数对CST方法的拟合精度影响较大。适当提高Bernstein多项式阶数可有效提高CST方法的拟合精度，但高阶Bernstein多项式容易导致拟合过程呈病态化。当W≥5时，参数化误差最大为9.7×10-4，已达到精度要求，因此本文选取W=5。采用5阶Bernstein时，翼型上、下翼面的参数化模型均有6个Bernstein多项式待定系数。以上翼面为例，A1、A2主要控制翼型前缘半径的大小，A3、A4主要控制翼型相对厚度和相对弯度的大小，A5、A6主要控制翼型后缘厚度的大小。通过CST方法对翼型几何特征进行参数化拟合，构建出计算机可识别的参数化模型，如图4所示。

图4 翼型参数化模型Fig.4 Parameterized model of airfoil

3.1.2 翼型气动特性敏感性分析

采用Xfoil软件求解翼型的气动特性。计算条件为稳态不可压缩气流，雷诺数Re=1.0×106，来流攻角θ=5°。

令各设计变量的均值为其标准差的10倍，各概率密度分布函数不完全相同，并采用NT-net法对全部设计变量分别进行抽样，构建各设计变量样本空间。总变量维数b=12，应用生成矢量产生GLP集合。为便于快速计算，根据前文所述Bernstein多项式待定系数的影响机制及翼型几何外形的平顺性要求，将变量维数分为4组，每组维数为3，查表得到生成矢量为（101；
1，40，85），则GLP集合为{gk=(gk1，gk2，gk3) 。}

采用改进的Morris法进行翼型气动特性的全局敏感性分析。以θ=5°时S832翼型的升阻比ε为输出量，设计变量zi（i=1，2，…，12）为输入量，其中z1，z2，…，z6和z7，z8，…，z12分别为依次对应的上、下翼面的设计变量。该12个设计变量是最小二乘法曲线拟合的数据点，依次控制5阶Bernstein多项式中的12个待定系数。以升阻比对设计变量的变异系数ζ为敏感性指标，采用改进的Morris法进行全局敏感性分析。θ=5°时升阻比对各设计变量的变异系数如图5所示。

图5 θ=5°时升阻比对各设计变量的变异系数Fig.5 Variation coefficient of lift-drag ratio to design variables at θ=5°

由图5可知，上、下翼面的升阻比对各设计变量的敏感性并不相同，其中对z9、z4、z1、z7、z6的敏感性较高。z4、z9控制待定系数A4、A9，间接控制翼型最大相对厚度和最大相对弯度的大小，表示翼型气动特性主要受到翼型最大厚度和最大弯度的影响，其中以z9的变异系数最大，表明下翼面翼型最大相对厚度和最大相对弯度对翼型气动特性的影响比上翼面翼型最大相对厚度和最大相对弯度对翼型气动特性的影响大；
z1、z7控制着待定系数A1、A7，间接控制翼型的前缘半径大小，表示翼型气动特性容易受到翼型前缘半径的影响；
z6控制待定系数A6，间接控制翼型后缘厚度大小，表示翼型后缘厚度也会影响翼型气动特性，但比前缘半径的影响程度要低。其余设计变量的变异系数均较小，表示升阻比对这些设计变量的敏感性较低，即这些设计变量对翼型气动特性的影响程度较低。

3.2 MCS验证

采用MCS法，每次抽取101个样本，共计12次，并建立翼型参数化模型，调用Xfoil软件求解翼型气动特性，然后采用Morris法进行全局敏感性分析。相对于抽样和建模时间，调用Xfoil软件及进行全局敏感性分析的时间可以忽略。在同等条件下，采用MCS法和NT-net法进行1次抽样的时间分别为205 s和180 s。采用MCS法进行全局敏感性分析所需的样本计算总数为n（i+2）［25］，其中n=101，i=12，则采用MCS法进行全局敏感性分析的计算次数为101×（12+2）=1 414，而NT-net法的计算次数为 101×12=1 212。MCS法和NT-net法的建模时长及进行全局敏感性分析的计算次数如表1所示。

表1 MCS法和NT-net法的建模时长及进行全局敏感性分析的计算次数Table 1 Modeling time of MCS method and NT-net method and calculation times of global sensitivity analysis

由表1可知，相比于MCS法，采用NT-net法的建模时间减少了300 s，总计算次数减少了202次。因此，采用NT-net法可以有效提高翼型气动特性敏感性分析的效率。

3.3 主要设计变量对翼型外形及气动特性的影响

在来流攻角θ=5°的工况下，翼型的气动特性主要受翼型几何特征的影响。在基于NT-net法得到的各设计变量样本空间中，分析主要设计变量对翼型几何外形及气动特性的影响。当设计变量变化时，翼型的最大相对厚度和最大相对弯度几乎同时发生变化，因此本文不对翼型最大相对厚度和最大相对弯度进行单独分析。后缘厚度对翼型气动特性的影响比前缘半径的影响小。因此，本文仅对前缘半径、最大相对厚度和最大相对弯度对翼型外形及气动特性的影响作深入分析。

3.3.1 前缘半径对翼型外形及气动特性的影响

选取样本空间中间接控制上、下翼面前缘半径大小的设计变量的最大值z1max、最小值z7min得到对应的待定系数A1、A7，其余待定系数与翼型S832一致，拟合出光滑的几何外形，命名为翼型1，如图6所示。

图6 以z1max、z7min间接拟合的翼型1的几何外形Fig.6 Geometric profile of airfoil 1 indirectly fitted with z1max、z7min

由图 6可知：当z1、z7分别达到z1max、z7min时，翼型上、下翼面对应数据点间的距离比S832的大，所控制的待定系数的变化幅度也比翼型S832的大，因此翼型1的前缘半径比翼型S832的大；
针对前缘区域，上翼面的变化幅度比下翼面的变化幅度大；
z1、z7对翼型其余部位的影响程度较小。

将翼型1的几何外形坐标导入翼型设计软件Profili和Xfoil，得到翼型1的前缘半径和升阻比。其与翼型S832的对比如表2所示。

表2 翼型1与翼型S832前缘半径和升阻比的对比Table 2 Comparison of leading edge radius and lift-drag ratio between airfoil 1 and airfoil S832

由表2可知，翼型1的前缘半径比S832增加了0.05 mm，升阻比提高了0.659 1。因此，当前缘半径适当增大时，翼型的升阻比提高，气动特性增强。因此，设计时翼型的前缘半径可适当取较大的值。

3.3.2 最大相对厚度和最大相对弯度对翼型外形及气动特性的影响

选取样本空间中控制上、下翼面最大相对厚度和最大相对弯度的设计变量最大值z4max、最小值z9min得到对应的待定系数A4、A9，其余待定系数与翼型S832的一致，拟合出光滑的几何外形，命名为翼型2，结果如图7所示。

图7 以z4max、z9min间接拟合的翼型2的几何外形Fig.7 Geometric profile of airfoil 2 indirectly fitted with z4max、z9min

由图7可知：当z4、z9分别达到为z4max、z9min时，翼型2上、下翼面对应数据点间的距离比翼型S832的大，所控制的待定系数的变化幅度也比翼型S832的大，因此翼型2的最大相对厚度和最大相对弯度比翼型S832的大；
针对翼型的中部区域，下翼面最大相对厚度和最大相对弯度的变化幅度比上翼面的大；
z4、z9对翼型其余部位的影响较小。

同理，可得到翼型2的最大相对厚度、最大相对弯度和升阻比。其与翼型S832的对比如表3所示。

表3 翼型2与翼型S832最大相对厚度、最大相对弯度和升阻比的对比Table 3 Comparison of maximum relative thickness,maximum relative camber and lift-drag ratio between airfoil 2 and airfoil S832

由表3可知，翼型2的最大相对厚度和最大相对弯度分别比S832增加了1.42%、2.47%，升阻比减小了1.096 4。因此，当最大相对厚度和最大相对弯度增大时，翼型的升阻比减小，气动特性减弱，这是由于来流攻角较小的缘故。因此，在来流攻角较小的工况下，翼型最大相对厚度和最大相对弯度可取较小值。

1）分析了Bernstein多项式阶数对CST参数化翼型设计中拟合精度的影响。在参数拟合过程中，适当增大Bernstein多项式的阶数可有效提高CST方法的拟合精度。当阶数W=5时，采用CST方法可以达到拟合精度的要求。

2）基于NT-net法构建了翼型参数化模型。阐述了NT-net法的计算原理，基于NT-net法抽样并构建了翼型参数化模型的设计空间。经MCS法验证，采用NT-net法可以有效提高翼型参数化表达过程中不确定性传播的拟合精度及效率。

3）通过NT-net法与Morris法的结合，确定了影响翼型气动特性的主要几何因素并提出改进策略。影响翼型气动特性的主要几何因素依次为最大相对厚度和最大相对弯度、前缘半径和后缘厚度等，其中当来流攻角较小时，前缘半径适当取较大值，最大相对厚度和最大相对弯度取较小值，可以有效提高该工况下翼型的气动特性。