参数化艇身阻力特性的全局敏感度及设计空间

2021-07-01王晓亮

哈尔滨工业大学学报 2021年6期

张宇，王晓亮

(上海交通大学航空航天学院，上海 200240)

平流层一般指海拔高度处于10～55 km的空间，由于平流层的垂直对流效应小，气流稳定，适合布置平流层飞行器完成诸如中继通信、监察和运输等任务，具有极大的军事及民用价值[1].平流层飞艇是一种能定点飞行、效费比高的平流层飞行器，其结构一般包括艇身、吊舱、尾翼和支撑骨架(如图1所示).对平流层飞艇而言，是否具备长航时是一项重要性能指标[2-3].平流层飞艇飞行过程中受到的阻力由动力系统平衡，阻力系数与动力消耗近似为正比关系，加之飞艇表面积巨大，驻留风阻很高，因此即便阻力系数略微减小，也能极大的减小动力系统的能量消耗，从而提高执行任务时间.通常，艇身阻力占到了全艇阻力的约2/3[4-5]，因此找到使艇身具有最小阻力系数的外形具有十分重要的意义[6-7].

图1 飞艇结构示意

目前，飞艇减阻一般采用优化方法进行. Lutz等[8]通过在艇身布置周向分布的点源/汇从而得到了压力场和速度场，并利用边界层模型得到了不同雷诺数下的最小阻力外形.Wang等[9]采用势流-边界层耦合方法和混合遗传算法对平流层飞艇艇身外形进行了优化. Geruti等[10]基于粒子群优化算法(PSO)提出了适用于考虑附加质量的非常规布局飞艇的优化框架. Kanikdale等[11]采用GNVR作为飞艇基础外形，提出了飞艇外形的多变量约束方法，并用模拟退火算法对外形进行了优化. Wang等[12]将气动、结构、能源和质量作为优化对象，通过综合目标方程对飞艇外形进行了多学科优化.一般地，飞艇艇身外形廓线可由八参数、四参数、GNVR、椭球形和“海豚”形描述[13-15].但这些方法能包含的外形范围较小，表达过于数学化，效率较低. 翼型参数化方法较上述飞艇外形描述方法丰富很多，能表征更为多样的外形空间.翼型参数化方法主要有Hick-Henne Bump Functions、B-splines、PARSEC和Class/Shape function Transformation (CST)等4种.其中，Hick-Henne Bump Functions方法通过叠加形状函数并改变形状函数因子αi对已有的基础外形进行扰动，例如: Zhong等[16]以RAE2822为基础翼型，通过引入10个形状函数来表示新的翼型；Yang等[17]利用开源软件SU2，引入5个形状函数分别对NACA0012和RAE2822翼型进行了优化；总的来说，Hick-Henne Bump Functions方法适用于对已有外形进行优化，难以形成设计空间.B-Splines方法通过在翼型周围改变控制点位置来更新翼型形状，例如Pérez等[18]使用B-spline方法控制螺旋桨不同展向处叶素的外形，通过设定不同的控制点数和容错系数得到不同的三维螺旋桨CAD模型；Martin等[19]应用扩展的B-spline方法对三维机翼形状进行了优化；该方法具有很强的凸包性和稳定性，但上述控制点没有具体的物理含义，不便于形成明确的设计参考. CST方法通过类函数C(x)和形状函数S(x)的乘积定义外形曲线，例如Masdari等[20]应用CST方法描述翼型形状，结合离散涡方法，以功率系数为目标，对Savonius涡轮进行了优化. 虽然该方法的适应性较强，但表达过于数学化，同样不利于形成明确的设计参考. PARSEC方法通过11个参数控制翼型形状，Chen等[21]通过PARSEC方法描述翼型形状，以功率系数为目标对垂直轴流风力机进行了优化，使功率系数比NACA0015翼型的结果提高了13.26%；PARSEC方法每个参数均有清晰的物理含义，利于形成明确的设计空间. Sripawadkul等[22]系统地研究了上述4种翼型参数化方法的特性，并给出了每种方法在不同指标下的评价结果，其中的正交性指标体现了各参数化方法中参数与外形是否一一对应，该指标在形成设计空间的过程中十分重要，结果表明只有PARSEC方法与CST方法能完全满足正交性.考虑到CST方法表达过于繁复，可选PARSEC方法来描述飞艇外形.

一般而言，优化过程是面向所有参数的，这往往导致在次要参数上会消耗许多不必要的时间.因此提取对艇身阻力系数最敏感的参数，降低优化设计维度，并由这些参数形成飞艇艇身外形设计空间，从而给予艇身初期设计一定的参考是十分必要的.基于此，本文以艇身阻力系数为目标，通过PARSEC参数化方法、CFD方法和Sobol全局敏感度分析方法得到对阻力系数最敏感的参数，进而形成了由这些参数构成的飞艇艇身外形设计空间.

1 数值方法及其验证

1.1 数值方法

艇身阻力系数可由势流理论、风洞实验和计算流体力学(CFD)方法得到[23].风洞实验的高耗时和高成本不适用于本文工作，针对飞艇发展的势流理论通常只对细长体预测的较准确，且黏性作用往往使势流理论的结果不准确.随着计算机技术的发展，CFD成为解决气动优化设计和流固耦合问题的重要手段.低速不可压缩流动的NS方程为

(1)

式中：ρ为流体密度;p为压力;μ为动力黏度;u为速度矢量.对式(1)进行无量纲化可得

(2)

选取Spalart-Allmaras模型求解湍流流动，Spalart-Allmaras模型适合求解航空外流场问题，其计算开销低，稳定性好[24].求解器采取不可压缩形式的压力基求解器.压力-速度耦合格式为“coupled”，空间梯度离散方法为“least squares cell based”，压力项采用二阶格式离散，动量和修正湍流黏度采用二阶迎风格式离散.Spalart-Allmaras模型通过求解关于修正湍流黏度的输运方程获得流场信息:

(3)

1.2 网格无关性验证

因艇身外形一般为旋成体，故简化为二维轴对称模型.艇身长度为L，艇身与入口和出口的距离均为60L，远场边界与旋转轴的距离也为60L，如图2所示.模型经过归一化处理，模型长度L取为1.边界条件设置见表1.流场的数值求解精度主要由y+和艇身沿轴向的网格种子数量决定.设置艇身表面的第1层网格高度为10-6L，计算得到的流线型旋成体模型轮廓线上的y+分布如图3所示，横坐标为归一化的艇身轴向位置，可见本文设置能使y+满足要求.如图4所示，对比了在不同Rev下艇身轴向网格种子数量与体积阻力系数Cdv的关系，可见当艇身轴向网格种子数量取为600时，Cdv几乎不再变化.后续数值分析均基于以上述网格划分模式进行.

图3 不同Rev下Model 4154表面的y+分布

图4 网格无关性检验

表1 边界条件类型

图2 计算域几何示意图

文中以飞艇体积的立方根值作为特征长度，因此Rev定义如下：

(4)

式中，V为飞艇体积.

总阻力系数Cdv由压差阻力系数Cdpv和黏性阻力系数Cdfv构成：

(5)

Fd=Fp+Ff,

(6)

Cdv=Cdpv+Cdfv.

(7)

式中：Fd为总阻力;Fp为压差阻力;Ff为黏性阻力.

1.3 算法验证

文献[25]中对于流线型旋成体做了详细的流动实验分析记录，在此选取6组模型作为验证对象，模型材质为红木，长度统一为2.74 m，通过在水洞中拖曳模型获取阻力，模型中轴线距离水面深度为2.74 m，模型标签名分别为4 154、4 156、4 158、4 162、4 164和4 175，对应长细比分别为4、6、8、7、7和5.流体介质的密度和动力黏度分别为998.2 kg/m3和0.001 003 kg/m·s，确保在数值计算中使CFD与实验的雷诺数相等. 对比结果如图5所示，计算和实验结果的平均相对误差为1.5%，最大相对误差为3.8%，满足工程误差许可，表明本文数值方法是可靠的.

图5 数值方法验证结果

为支撑数值方法中雷诺数能表征流动情况这一结论，现从具体算例上进行结果分析.以4 154模型为对象，取流动雷诺数Rev为5×106，固定来流速度U∞为10 m/s，通过调整流体密度ρ与流体动力黏度μ使Rev保持不变.各种参数配合下4 154模型受到的阻力情况见表2，可见当雷诺数Rev一定时，模型的阻力系数可视为常数.因此不论从控制方程还是计算方法来看，雷诺数均能表征流动情况.

表2 Rev=5×106时的阻力系数

2 飞艇外形参数化方法

PARSEC方法通过11个参数控制翼型形状，每个参数均有清晰的物理含义.鉴于飞艇艇身的旋成体特性，只需选取上半部分的参数，于是可通过8个参数描述飞艇外形，如图6所示.

图6 基于“PARSEC”的8参数飞艇艇身外形

飞艇外形可由6阶多项式表达：

(8)

式中的待定系数可通过求解下式获得[21]，xte为弦长：

(9)

3 全局敏感度分析方法

敏感度体现了变量自身的改变对系统的影响程度.一般地，敏感度分析方法分为局部敏感度分析和全局敏感度分析.例如，龙卷风图,基于分化的方法和筛查法属于前者;基于回归的方法,基于方差的方法,转换不变方法和基于密度的方法属于后者[26].局部敏感度能体现输入空间一点附近的特性，全局敏感度还考虑了参数之间的相互影响，结果更合理可靠.

3.1 Sobol全局敏感度分析方法

假设物理模型能被方程f(x)描述，输入x=(x1,x2,…,xn) 是n维空间的一点且xi(i=1,2,…，n)遵循[0,1]上的均匀分布.全体x构成一个n维立方体:

Cn={x|0≤xi≤1;i=1,2,…,n}.

(10)

输出项f(x)能被分解为

f1,2,…,n(x1,x2,…,xn),

(11)

式中f0为常数，且f1,2,…,s(x1,x2,…,xs)关于它自身变量的积分为0.

(12)

积分式(11)有

(13)

由于式(11)右侧至少有一个变量是不同的，故具有正交性，再由式(12)可知:

0,(i1,i2,…,is)≠(j1,j2,…,jl).

(14)

对式(11)积分和平方,即

(15)

式(15)右侧第2项被称为总方差，写为

(16)

偏方差定义如下:

(17)

偏方差和总方差的关系为

(18)

引入比率S1,2,…,s,即

(19)

对于具有n个输入变量的离散样本，本文需要计算每个S1,2…,s的值，但对于需要借助CFD软件才能获得的输出，这样势必造成极大的时间开销.基于此，本文将所有输入变量分为两个子集，子集Xi仅仅包含变量xi, 另一子集Xei包含除了xi的所有变量. 因此总方差可写为

D=Di+Dei+Di,ei,

(20)

引入新的参量STi，即

STi=Si+Si,ei=1-Sei,

(21)

Si表征了变量xi单独对总方差的影响，STi体现了变量xi对系统的“总体”影响，即不仅仅包含变量xi自身也包含了变量xi和余下变量的任意组合对输出的影响.因此，本文称Si和STi为变量xi的一阶敏感度和全局敏感度.

3.2 离散数据的敏感度分析

事实上，本文经常面对的是离散的输入与输出数据，其中没有明确的函数关系.现在假设有一具有n个输入变量的系统，抽样产生N个样本.计算之前，本文需要准备两组输入数据x1,x2,…,xn和x′1,x′2,…,x′n，写成矩阵形式如下：

(22)

有关量可由蒙特卡洛方法近似得到[27]：

(23)

(24)

(25)

(26)

式中：fi,j为将变量xj替换为x′j所得到的输出；f′i,j为将变量x′j替换为xj所得到的输出.于是本文能通过上述算法和定义获得输入变量的一阶或全局敏感度.

4 敏感度分析和艇身阻力特性研究

4.1 敏感度分析

表3给出了基于PARSEC方法的飞艇外形参数范围.头部半径rh不超过艇身长度1；最大半径rd处于0.05～0.25之间，对应艇身长细比处于2～10之间，能涵盖大多数常规艇身；最大半径位置xd处于0.3～0.5之间，以避免艇身最大截面位置过于靠前或靠后，以此减少艇身外形出现畸变的概率；尾部张开角βt处于20°～40°之间，以保证艇身尾部光滑收缩；尾部偏移角αt，尾部厚度Δyt和尾部高度yt被设定为0，以确保表征的外形为旋成体.如图7所示，不恰当的参数组合会导致病态外形出现，实线外形包含了负值，虚线外形有多个极值，因此这些病态外形要在抽样时被剔除.基于表3中的参数范围，得到所有非病态外形样本构成的范围如图8所示，可见所取的参数能很好地保证艇身廓线的多样性.

图7 不合理的参数组合造成的病态外形

图8 艇身廓线构成的样本空间

表3 飞艇外形参数范围

飞艇飞行高度设定为20 km，当地空气密度为0.088 9 kg/m3，空气动力黏度为1.421 61×10-5kg/m·s，大气压强为5 529.31 Pa.一般地，平流层飞艇的绕流雷诺数在6.0×106和1.2×107之间，对一典型体积为3.0×105m3的飞艇而言，当来流速度为20 m/s时，其雷诺数为8.37×106.故本文取雷诺数为8.0×106进行后续研究，进行数值分析时调整流体密度使所有工况的绕流雷诺数一致.

图9显示了每个外形参数对各阻力系数的总敏感度.可以发现，参数rh和rd对总阻力系数的敏感度分别达到68.8%和66.0%，除了xd以外的其余参数对总阻力系数几乎没有影响，且xd的总敏感度仅为2.1%，故剩余参数可设为确定值.观察由参数rh和rd组成的设计空间，该空间可近似看作被3条线段和1条三次曲线包围形成，如图10所示.该空间的数学表述为

图9 不同参数的总敏感度

图10 病态和设计空间

(27)

4.2 阻力特性分析

分别设置参数xd，y″d和βt的值为0.35，-0.5和30°.然后将设计空间离散，由同样的计算条件得到参数rh和rd与各阻力系数的关系，如图11所示，随着设计空间向右上方推进压差阻力系数逐渐升高，而黏性阻力系数与此趋势恰好相反，二者的综合效果使总阻力系数在设计空间上存在极小区域.由此可得到参数rh和rd的最佳组合，另外由图12可知，在其他参数确定的情况下随着xd的增加，总阻力系数先减小再增大，当xd=0.435时，总阻力系数取最小值.

图11 rh和rd形成的关于各个阻力系数的云图

图12 Cdv和xd的关系

综上所述，可获得使艇身具有最小阻力系数的参数组合，该艇身外形的头部半径rh为0.012 2，最大半径rd为0.95，最大半径位置xd为0.435.上述所得艇身长细比为5.263，与文献[28]中的结论吻合.在实际设计阶段，头部半径rh的范围可取为[0,0.06]，最大半径rd的范围可取为[0.08,0.12]，该范围内的阻力系数大小浮动程度较小，均可较好的降低艇身阻力系数.