基于POD和代理模型的热气防冰性能预测方法

2023-01-31杨倩郭晓峰李芹董威

航空学报 2023年1期

杨倩，郭晓峰，李芹，董威

上海交通大学机械与动力工程学院，上海 200240

在结冰气象条件下飞行时，飞机部件表面会产生积冰，严重危害飞行安全［1-2］。近些年，飞机结冰及由此引发的飞行事故越来越受到关注［3-4］。民用航空飞行器广泛采用热防冰系统将部件表面加热至0 ℃以上，以达到防冰目的。其中，热气防冰系统是大型民用飞机上常见的热防冰手段。从压气机某级引出的热空气经由管道输送到机翼、进气道等部件前缘，再经由笛形管表面开有的射流孔流入防冰腔，并沿蒙皮内表面流动，将热量通过蒙皮传导到外表面，保证机翼、进气道等部位不发生结冰。

防冰腔内部结构是影响防冰表面温度、溢流水分布等性能的主要因素。在防冰系统优化设计过程中需对内部结构进行多轮迭代，以保证在一定引气温度和热气流量下达到最优的防冰效果和最高的能量利用效率。评估防冰系统性能的方法有冰风洞试验、数值仿真等。其中冰风洞试验准备周期长，实施成本高，缝翼、全尺寸短舱等部件还需进行模型缩比或试验条件缩比［5-6］，在设计阶段不适用。利用数值仿真评估防冰性能包括蒙皮内外流场计算［7-8］，防冰性能计算［9］等多个过程，迭代求解速度仍旧较慢。特别是当采用遗传算法等方法开展优化设计时，需要上千、甚至上万次性能评估，数值仿真所需的时间成本也是不能接受的。

代理模型是基于给定的描述客观事实的数据集进行插值或者回归构造出的近似替代模型，用以表达复杂系统，被广泛应用于多个研究领域［10-14］。代理模型需要比原始模型计算量小、计算时间短、且能保证计算精度不显著降低。代理模型的结构大多是多变量输入，单目标值输出。例如输入防冰腔内笛形管射流孔间距、射流孔角度，输出防冰表面平均温度或总溢流水流量。但是，在热气防冰系统优化设计过程中仅仅知道表面平均温度或总溢流水流量并不能有效指导结构参数优化方向。全面描述防冰系统性能需要建立的代理模型个数为防冰表面总网格点数和选定防冰性能参数总个数的乘积，这无疑会增加构建代理模型的成本。

本征正交分解（Proper Orthogonal Decom⁃position，POD）能够将复杂系统分解为若干个基模态来捕捉系统的主要特征，并利用线性叠加基模态实现复杂系统重构［15-18］。利用本征正交分解能够将三维防冰表面温度、溢流水等性能参数分解为表达其特征的基模态，并得到利用基模态重构防冰性能所需的线性拟合系数，从而将复杂的三维表面性能参数样本降阶为可以用于代理模型训练的拟合系数样本。但是，防冰表面溢流水分布受到外部结冰环境、内部防冰结构的多重影响，流动存在干、湿间断分区，需要合理选择基模态阶数、拟合系数维度才能实现高效的性能预测。

基于本征正交分解和代理模型提出了一种实现热气防冰性能快速预测的方法。首先在热气防冰结构设计变量空间内进行拉丁超立方抽样［19］（Latin Hypercube Sampling，LHS），得到均匀分层的笛形管结构参数样本，并求解结构参数样本对应表面温度、溢流水分布等防冰性能。之后，利用本征正交分解找出表达防冰性能主导特征的基模态以及重构防冰性能所需拟合系数。最后，基于支持向量机回归［20-21］（Support Vector Regression，SVR）机器学习方法建立结构参数样本和拟合系数间的代理模型，实现笛形管热气防冰优化设计过程中对防冰性能的预测与评估。针对三维缝翼笛形管热气防冰系统所开展的验证表明该方法对热气防冰表面温度和溢流水均有较好的预测效果。

1 本征正交分解方法

本征正交分解可以在观测所得复杂物理场的基础上获得一系列基模态。借助Sirovich［22-24］引入的方法，令线性无关向量集合中的每一个元素均为l维空间Ω∈Rl中的一个向量，称之为“快照（Snapshot）”。找出快照所张空间Ψ中的一组规范正交基，使得集合中的元素在这组基上的投影最大：

式中：( · ， · )表示内积运算。规范正交基中的基向量可以用向量之间的线性叠加来表示，即求得系数即可求得规范正交基Φ(j)。对于规范正交基Φ(j)，其第j个基模态中的第i个系数为快照间协方差矩阵Rn×n的第j个特征值对应特征向量的第i个元素，即：令λ(1)，λ(2)，⋅⋅⋅，λ(n)为矩阵Rn×n的从大到小排列的n个互异特征值，其中Rn×n的第i行第j列元素为为λ(i)对应的特征向量。第i个POD基模态可以表示为

n个快照可以分解得到n个基模态。理论证明，保留少量POD基模态就可以有效捕捉原始数据的主要特征。特征值的大小表征了POD基模态对原始数据特征的代表程度，前m个基模态所包含的广义能占总能量的比率为

可以按照特征值λ由大到小衰减、能量占比E由小到大增加的序列选择适当的POD基模态阶数。

2 支持向量机回归模型

支持向量机是基于统计学理论的机器学习方法，可以用于回归分析。该方法使用结构最小化原则，可以更好地解决小样本学习问题，同时使用核函数将非线性问题转化为线性问题，降低了算法的复杂程度。假定训练样本集首先考虑用线性回归函数：

式中：w和b为回归函数f(x)的权值和阈值。采用最小化欧几里德空间的范数来寻找最小的w。假定所有训练数据(xi，yi)都可以在精度ε下用线性函数拟合，那么寻找最小w的问题就可以表示成凸优化问题：

约束条件：

考虑允许存在拟合误差，引入松弛因子ξi≥0和将优化问题转化为

约束条件：

式中：i=1，2，…，n；常数C>0是预先给定的，用来平衡回归函数f的平坦程度和偏差大于ε的样本点个数，函数的图形如图1所示。

图1 支持向量机回归Fig. 1 Support vector regression

将上述优化问题转化为相应的对偶问题，同时引进核方法将对偶问题转化为求解以下约束问题的最大值，解得Lagrange乘子

式中：K( · ， · )表示核函数。约束条件为

出于稳定性考虑，b的求解采用支持向量的平均值，其中

得到目标回归方程：

3 基于POD和支持向量机回归的防冰性能快速预测方法

基于上述本征正交分解和支持向量机回归模型提出了一种全新的热气防冰性能快速预测方法。首先在笛形管设计变量空间进行拉丁超立方抽样，得到n个q维的设计参数向量P(i)=[xi1，xi2，…，xiq]，i=1，2，…，n。针对P(i)依次开展防冰性能数值仿真，得到n个防冰性能样本。将固体蒙皮外表面温度和溢流水分布结果整理成快照的形式：

快照集合U的第i列表示输入结构参数P(i)对应的l个网格点上的温度值（T1i~Tli）和溢流水流量值

采用奇异值分解求解快照集合U的基模态及基模态对应特征值。对U进行奇异值分解可以得到

式中：W为一个2l×2l的矩阵；Σ为一个2l×n的矩阵，其主对角线上元素称为奇异值，主对角线以外的元素全为0；V为一个n×n的矩阵。W和V都是酉矩阵，满足WTW=I，VTV=I。W中的每一列都是快照集合U的一个基模态，Σ中主对角线上每一个元素的平方都是对应基模态的特征值λ。截取前m个基模态就可以捕捉到U的绝大部分特征。通过线性叠加这m个基模态可以实现对U的近似拟合：

其中：U(i)表示快照集合U的第i列，且m

式中：M为m×m的对称矩阵为m维向量由于W是酉矩阵，且为W矩阵的列向量，满足Mij=代入式（17）得至此，输入结构参数向量P(i)所确定的防冰性能样本就可以用m维系数向量来实现近似拟合。

采用代理模型建立输入参数P(i)和拟合系数b(i)之间的关系就可以预测已给出训练样本之外的输入参数对应的防冰性能。当截取前m个基模态时，需要针对m维的拟合系数建立m个代理模型，即

完成代理模型训练后，依次向m个代理模型输入测试样本向量1，2，…，ntest，预测得到测试样本对应的m维拟合系数b(k)，如图2所示。将拟合系数向量b(k)和截取的基模态输入式（16），就可以实现P()k对应防冰性能的预测，即防冰表面l个网格点上的温度分布和溢流水流量分布。

图2 代理模型的输入和输出Fig. 2 Inputs and outputs of surrogate models

4 热气防冰表面温度和溢流水快照构建计算方法

传统防冰性能数值仿真方法求解热气防冰表面温度、溢流水分布等性能快照时，同时计算蒙皮外部冷空气流场和内部热气射流流场，工作量非常大。为了加快性能快照积累速度，将冷空气来流和过冷水滴撞击简化为蒙皮外表面边界条件，将由给定笛形管结构参数决定的热气射流流动简化为蒙皮内表面边界条件。仅针对蒙皮固体计算域，在其外表面建立溢流水质量和能量守恒方程，迭代计算防冰热载荷至表面温度和溢流水均收敛。

4.1 蒙皮外表面防冰热载荷计算

蒙皮外表面空气流场计算模型与选取研究模型保持一致，将模型位于流场计算域中间。近壁面网格进行局部加密处理，以满足Spalart-Allmaras湍流模型对边界层的计算要求。空气流场计算域入口和出口分别为速度入口和压力出口边界条件。

求解防冰部件表面局部水滴收集系数是开展防冰热载荷计算的前提。采用欧拉法求解水滴相场［25］，水滴相控制方程为

式中：α表示当前控制容积内水滴相体积分数；ρa、ρd分别表示空气、水滴的密度；ua表示空气运动速度；ud表示水滴运动速度；d表示水滴平均直径；CD为水滴阻力系数。

当水滴相控制方程计算收敛时，部件表面的水滴体积分数能够反映当地液态水含量的实际情况，从而计算得到当地局部水滴收集系数β：

式中：u∞表示未受到扰动位置的水滴速度；α∞表示未受到扰动位置的水滴相体积分数；n表示撞击位置的单位法向量。

防冰系统工作时撞击的过冷水滴升温并蒸发，没有蒸发的随着气流方向溢流形成水膜，防冰表面溢流水控制容积中的各质量项和能量项如图3所示。

图3 防冰表面控制容积质量与能量守恒示意图Fig. 3 Schematic diagram of mass and energy conser⁃vation of control volume on anti-icing surface

考虑图3所示防冰表面溢流水控制容积，建立质量守恒方程［26-27］：

根据防冰表面能量平衡，建立图3所示控制容积能量守恒方程［26-27］：

防冰系统工作时防冰热载荷与防冰系统提供的热量达到平衡，根据式（23）可以得到当前温度条件下蒙皮外部防冰热载荷Q̇load：

4.2 热气射流冲击换热计算

利用热气射流换热特性关联式来代替防冰腔内部流场数值仿真，从而获取给定笛形管结构参数样本所对应蒙皮内表面边界条件。根据热气射流的流动特性，可以将高温气体冲击蒙皮内表面划分为3个流动区域［8］，如图4所示。根据射流孔所在截面笛形管圆心、射流孔圆心、射流驻点三点共线在蒙皮内表面寻找射流驻点，并以驻点为基准沿径向求解蒙皮内表面的对流换热系数分布。

图4 热气射流流动示意图Fig. 4 Schematic diagram of hot air jet flow

适合的关联式是准确预测热气射流换热特性的关键。综合考虑笛形管表面射流孔孔径D、射流孔到蒙皮内表面的距离Zn、射流孔排布方式；热气射流出口雷诺数为单个射流孔热气质量流量，μa为气流黏性系数；不考虑热气射流之间相互作用［28］，选择由Goldstein等［29］利用试验整理得到的单孔射流冲击平均努塞尔数的计算关联式：

式中：R表示蒙皮内表面任意位置到射流驻点的距离；Nu表示当地努塞尔数；试验参数范围为61 000

在式（25）两侧对R求导，并整理得到

利用当地Nu数计算得到蒙皮内表面对流换热系数为

式中：h表示当地对流换热系数；ka表示气流热导率。

5 算例验证及分析

5.1 验证算例

为了验证提出的基于本征正交分解和支持向量机回归模型的热气防冰性能预测方法，选用某三维缝翼笛形管热气防冰模型开展验证计算。如图5（a）所示，模型由缝翼和翼盒2部分组成，弦长为1.2 m，展长为0.25 m，蒙皮厚度为2 mm。笛形管表面射流孔采用钻石型排布。如图5（b）所示，验证算例输入结构变量包括：① 展向2排射流孔间距L；② 笛形管轴线x方向坐标xpic；③ 笛形管轴线y方向坐标ypic；④ 中间排射流孔出流方向和水平方向夹角θ0；⑤ 两侧射流孔和中间排射流孔夹角θ1。因此，样本参数向量维度q=5，所有结构变量取值区间如表1所示。笛形管外径为60 mm，射流孔孔径为2 mm，均为固定值。由拉丁超立方抽样共生成均匀分层的1 000个训练样本和200个测试样本，部分训练样本和测试样本的结构参数分别如表2和表3所示。

图5 预测方法验证模型Fig. 5 Validation model for estimation method

表1 热气防冰结构参数取值区间Table 1 Range of hot air anti-icing geometric parameters

表2 训练样本Table 2 Training samples

在给定飞行条件和结冰条件下开展蒙皮外表面压力、气流速度、气流与表面对流换热系数、局部水滴收集系数计算，认为蒙皮外表面冷空气流场和水滴场不受内部笛形管结构参数变化的影响，仅需计算一次。空气速度为80 m/s，攻角为2.5°，环境温度为264.15 K，平均水滴直径为20 μm，液态水含量为0.5 g/m3。

利用热气射流冲击换热计算得到防冰结构参数对应蒙皮内表面对流换热系数，计算中输入笛形管结构参数：L、xpic、ypic、θ0、θ1；并给定热气射流总温459.15 K，单个射流孔流量与孔间距L有关，为1.33⋅(L/30) g/s。

表3 测试样本Table 3 Testing samples

在给定固体蒙皮内外表面边界条件的基础上，迭代计算蒙皮外表面温度和溢流水分布，直到两者均收敛，获得该结构参数样本对应防冰性能快照。用于开展防冰性能计算的蒙皮固体域网格为结构网格，共230 400个，如图6所示。蒙皮外表面共23 040个网格，即式（14）中l=23 040。

图6 缝翼蒙皮结构化网格Fig. 6 Structured mesh of slat skin

为了有效评价POD拟合结果、代理模型预测结果的优劣，定义温度、溢流水流量的拟合误差或预测误差为

式中：ERT、ERṁout为所有样本的平均拟合误差或平均预测误差；MAET，i、MAEṁout，i分别为第i个样本温度结果和溢流水流量结果的拟合误差或预测误差，具体定义为

现有飞机防冰计算中针对溢流现象采用的物理模型主要集中于溢流水控制容积质量与能量守恒分析，对溢流水膜破裂、演化规律及流动形态预测等方面还需要开展深入研究，来进一步提高数值仿真积累快照样本的精度。在针对建立的预测方法进行效果评价时，为了不引起混淆，将结构参数P(i)或P(k)通过数值仿真求解得到的结果称为“实际结果”，利用POD拟合得到的结果称为“拟合结果”，利用代理模型预测得到的结果称为“预测结果”。

5.2 本征正交分解拟合结果

基于本征正交分解可以得到防冰性能快照集合U的基模态以及基模态对应特征值。如图7所示，基模态对应特征值衰减迅速。根据式（4）计算得到了前m个基模态所包含广义能的比率，发现第1阶POD模态就包含了99.887%的能量，前10阶包含了99.998%以上的能量。

图7 前m个POD基模态对应特征值和广义能Fig. 7 Eigenvalues and energy captured by the firstmPOD basis

图8给出了利用不同阶数POD基模态得到温度、溢流水拟合结果和训练样本实际结果的平均绝对误差。可以看出，随着截取阶数增多，每一个训练样本拟合结果和实际结果间的平均绝对误差都大幅下降。但是对比发现，当截取阶数为20时，ERT已经<0.5 K，截取>80阶基模态时误差基本保持不变，而ERṁout只有在阶数>120时才基本不变，<0.005 g/（m·s）。

图8 拟合误差随基模态数量变化Fig. 8 Fitted error changing with number of basis truncated

图9 前m个POD基模态包含的广义能Fig. 9 Energy captured by the firstmPOD basis

图10 训练样本#1拟合结果与实际结果对比Fig. 10 Comparison between fitted results and true results of Training sample #1

依据POD基模态截取阶数对拟合结果和实际结果之间误差的影响，截取前160阶基模态重构防冰性能，部分训练样本重构结果如图10、图11所示。图中s表示从缝翼前缘测量的距离，正值表示缝翼上表面，负值表示缝翼下表面，s=0表示缝翼几何驻点；z表示缝翼展向方向。实线对应笛形管钻石型排布中双排孔位置，虚线对应单排孔位置，并对这两处位置的温度、溢流水流量的拟合结果（POD）和数值结果（CFD）进行了详细对比。对比防冰表面溢流水流量分布，拟合结果和实际结果吻合得也较好，不仅能够准确捕捉溢流水由水滴撞击区向缝翼后缘流动的趋势，而且能够准确预测溢流水流量数值。

图11 训练样本#2拟合结果与实际结果对比Fig.11 Comparison between fitted results and true results of Training sample #2

5.3 代理模型预测结果

综合考虑不同阶数基模态的拟合效果和代理模型训练所需成本，选择截取前160阶POD基模态用于防冰性能预测。因此，需针对建立160个代理模型。基于支持向量机回归方法训练代理模型的过程中，采用网格搜索（Grid Search）和交叉验证（Cross Validation）来获得最佳的核系数γ和惩罚参数C的组合。代理模型训练完成后依次预测测试样本P(k)对应的bk1，bk2，…，bk160，再利用式（16）线性叠加截取的160阶基模态，实现防冰性能预测。

图12展示了针对1 000个训练样本和200个测试样本采用代理模型预测得到的温度、溢流水预测结果与实际结果的平均绝对误差。对比发现，训练样本和测试样本误差量级相当，保证了所建立预测方法有较好的泛化性，可以用于未知笛形管结构参数对应防冰性能的快速预测。

图13~图15分别展示了代理模型预测结果与测试样本实际结果的对比。可以发现，预测温度结果和实际结果符合得较好，防冰表面各个位置的温度数值和变化规律都预测得较为准确。从溢流水预测结果与实际结果的对比来看，建立的预测方法在水滴撞击区域内能够实现较为准确的溢流水分布预测，如测试样本#1。但在水滴撞击区域之外与实际结果存在差异，如测试样本#2下表面和测试样本#3下表面。造成这一现象的原因可解释为：由于热气射流会形成局部高温区域，高温区域内液态水完全蒸发或大量蒸发，导致防冰表面出现沿展向的周期性干、湿分区或沿展向的周期性流量先增大后减小的溪流状水膜，即溢流水膜在该位置存在间断特征。虽然能够利用基模态捕捉到这类特征，但是在线性叠加过程中，基模态之间会将捕捉到的间断特征相互抵消，这一现象和翼型跨声速流场中激波预测存在误差相似［30］，只能通过选择合适的基模态阶数来平衡全场溢流水预测精度和间断处溢流水预测精度。

图12 训练样本和测试样本预测误差Fig. 12 Estimated error of training samples and testing samples

图13 测试样本#1预测结果与实际结果对比Fig. 13 Comparison between estimated results and true results of Testing sample #1

在Intel（R） Core（TM） i9主频3.70 GHz 的CPU上单线程开展200个防冰结构参数测试样本的性能预测共需要9.41 s，即单个测试样本仅需大约0.048 s，较需要数小时甚至数天才能完成的传统数值仿真计算方法有了大幅提升，这对采用遗传算法等优化方法开展热气防冰系统优化设计具有重要意义。