﻿ 基于降阶模型的翼型结冰冰形预测方法<sup>*</sup>
 文章快速检索 高级检索

Ice shape prediction method of aero-icing based on reduced order model
LIU Teng, LI Dong, HUANG Ranran, ZHANG Zhenhui
School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, China
Received: 2018-08-14; Accepted: 2018-12-28; Published online: 2019-01-16 15:14
Corresponding author. LI Dong.E-mail:ldgh@nwpu.edu.cn
Abstract: The numerical simulation prediction of aero-icing ice shape is usually complicated and time-consuming. In order to predict ice shape more quickly and accurately and to reduce computational resource consumption, a quick ice shape prediction method based on proper orthogonal decomposition and Kriging model is proposed in this paper. Using a database of high-fidelity CFD numerical simulations to build sample space, in view of the change of attack angle, the procedure for predicting the ice shape by reduced order model is introduced. With the data sampling method of parameter space, multiparameter prediction of the ice shape is achieved. Meanwhile, the related methods of the advanced Blind-Kriging model and the improvement of the prediction results are studied. The result shows that the prediction results of the airfoil icing shape using the reduced order model agree well with the CFD numerical simulation results. The conclusion is made that the reduced order model is a quick and accurate approach for predicting the ice shapes of airfoil icing.
Keywords: ice-shape prediction     reduced order model     proper orthogonal decomposition (POD)     Kriging model     Blind-Kriging model

1 基本方法

 图 1 降阶模型建立流程 Fig. 1 Reduced order model building process
1.1 参数空间采样

LHS方法是由Mckay等[12]在1979年提出的试验设计方法，该方法的基本原理是：如果需要得到n个样本点，则利用m个设计变量将参数空间均匀划分为等概率的n个区间，这样整个参数设计空间被分为nm个小区间，然后对其进行n次取样。设计空间抽样时一方面需要保证样本点在每个小区间上都是随机选取的；另一方面需要各个样本点在任一维度上进行投影，每个小区间中只会有一个样本[12]。通常取样的算法为

 (1)

 (2)

1.2 POD方法

POD方法是基于物理场原始数据快速准确地重构物理场，同时也可预测物理场在某些未知位置的数据的一种数值方法。POD方法的核心思想是利用特殊正交基函数的迭加来实现参数的近似[14]。它最初主要用于湍流拟序结构的分析研究中，POD如今已成功地应用于各类领域，如机器学习中的主成分分析(Principal Components Analysis, PCA)法，海洋学和气象学中的经验正交函数法，心理学和经济学中常用的要素分析法等[15]

POD的基本理论[16]为假设存在这样的物理场y=y(x, t), 在某个空间可以表示为如式(3)所示无穷级数:

 (3)

 (4)

 (5)

 (6)

 (7)

 (8)

 (9)

 (10)

1.3 Kriging模型

Kriging模型[17]是由南非工程师克里金(Krige)于20世纪50年代提出的, 现已广泛应用于各个领域，本文以普通克里金插值方法为例来说明克里金插值的原理[18]

 (11)

 (12)

 (13)

 (14)

 (15)

 (16)

β0σ2代入后，取对数忽略常数项，优化问题转化为

 (17)

2 冰形模拟及样本获取

 图 2 FENSAP-ICE软件模块相互关系 Fig. 2 Relationship among FENSAP-ICE software modules

 图 3 FENSAP-ICE冰形模拟结果 Fig. 3 FENSAP-ICE ice shape simulation results

 图 4 冰形模拟与实验对比验证 Fig. 4 Comparison and validation of ice shape simulation and experiment
3 降阶模型应用算例与结果分析 3.1 单参数冰形预测

 计算状态 数值 来流速度/(m·s-1) 100 液态水含量(LWC)/(g·m-3) 1 平均水滴直径(MVD)/μm 20 环境压力/Pa 101 325 结冰温度/K 263.15 结冰时间/s 360

 图 5 不同飞行迎角的冰形 Fig. 5 Ice shapes of different flying attack angles

 特征值序号 λx/10-6 λy/10-6 1 174 832 6 10 328.48 2 8.937 12 20.761 14 3 3.189 29 5.440 25 4 1.121 569 1.452 005 5 0.613 951 0.556 798

 图 6 单参数POD与CFD冰形对比 Fig. 6 Ice shape comparison between single-parameter POD and CFD
3.2 多参数冰形预测

 参数 最小值 最大值 飞行迎角/(°) 0 5 飞行速度/(m·s-1) 90 130 结冰温度/°F -22 +32 MVD/μm 15 40 高度/ft 0 22 000

 图 7 连续最大结冰条件[21] Fig. 7 Continuous maximum icing condition[21]

 图 8 样本空间分布 Fig. 8 Sample space distribution

LHS方法是一种“充满空间的设计”，选取的样本点可以均匀地散布于参数空间。

 算例状态 数值 飞行迎角/(°) 1.325 飞行速度/(m·s-1) 91.875 MVD/μm 25.5 高度/ft 12 031 结冰温度/°F 16 结冰时间/s 360

 图 9 多参数POD与CFD冰形对比 Fig. 9 Ice shape comparison of multiparameter' POD and CFD
4 插值模型 4.1 Blind-Kriging模型

 (18)

 (19)

Blind-Kriging最重要的步骤是识别出未知函数fi，这是通过函数(变量)选择技术在简单函数构成的候选函数集中选择得到，然后就可以得到预测模型的第一部分f(x)Tβm

X←samples

b1, …, bt{Candidate features}

b=Cte{Selected features}

θ0=maxθlikelihood(X, θ, b)

M0=construct(X, θ, b){Construct ordinary Kriging model}

α0=evaluateMeasure(M0){Assess accuracy}

i=0

while improvement(α){Accuracy improves?}

i=i+1

β=rank(b1, …, bt)

j=maxj(|βj|)

b=bbj

θi=maxθlikelihood(X, θ, b){Optional}

Mi=update(Mi-1, θi, b){Intermediate Kriging model}

αi=evaluateMeasure(Mi){Assess accuracy}

endwhile

θfinal=maxθlikelihood(X, θ)

Mfinal=update(Mi, θfinal, b){Final Blind Kriging model}

4.2 模型应用

 计算状态 数值 MVD/μm 20 环境压力/Pa 101 325 结冰温度/K 263.15 结冰时间/s 360

 图 10 冰形预测结果对比 Fig. 10 Comparison of ice shape prediction results
5 结论

1) 通过对POD预测结果与CFD结果进行对比可知：降阶模型可以得到较好的冰形预测效果，只是在冰角附近及结冰极限位置预测结果有细微的差别，表明该方法应用于翼型结冰冰形的预测是准确、可行的。

2) 利用贝叶斯框架下的Blind-Kriging模型进行数据插值虽然需要更多的时间完成模型建立，但可以改善降阶模型的精度，可以取得更为理想的冰形预测效果。

3) 虽然样本的取得需要一定的计算时间，但是降阶模型建立完成之后，便能够快速、准确地预测一定范围内任意状态的冰形。

4) 本文方法对涉及大参数研究的应用极具价值，如果预先计算足够的快照，则通过构建POD/Kriging逼近而不是使用全部的CFD解决方案来进行参数研究可以很容易节省大量时间。这为结冰试飞认证等工程应用提供了一种有效可行的方法，同时该方法可以推广应用至翼型设计阶段，对于容冰设计中的翼型结冰影响分析以及容冰优化设计是一种可靠准确的预测手段。

 [1] 蒋天俊.结冰对飞机飞行性能影响的研究[D].南京: 南京航空航天大学, 2008. JIANG T J.Investigation of icing accretion influences on aircraft flight performance[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2008(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=D052562 [2] 姚若鹏.翼型的结冰数值模拟及相关控制研究[D].南京: 南京航空航天大学, 2012. YAO R P.The numerical simulation of ice accretion on airfoil and control research[D].Nanjing: Nanjing University of Aeronautics and Astronautics, 2012(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=D280454 [3] 杨胜华.二维飞机结冰过程仿真[D].北京: 北京航空航天大学, 2010. YANG S H.Two-dimensional in-flight ice accretion simulation[D].Beijing: Beihang University, 2010(in Chinese). [4] 申晓斌, 郁嘉, 林贵平, 等. 基于特征正交分解法的翼型结冰冰形快速预测[J]. 航空动力学报, 2013, 28(4): 807-812. SHEN X B, YU J, LIN G P, et al. Fast prediction of ice shape based on proper orthogonal decomposition method[J]. Journal of Aerospace Power, 2013, 28(4): 807-812. (in Chinese) [5] GALLIVAN K, GRIMME E, VAN DOOREN P.Pade approximation of large-scale dynamics systems with lanczos methods[C]//Proceedings of the 33rd IEEE Conference on Decision and Control.Piscataway, NJ: IEEE Press, 1994: 443-448. [6] CAMUSSI R, GUJ G. Orthonormal wavelet decomposition of turbulent flows:Intermittency and coherent structures[J]. Journal of Fluid Mechanics, 1997, 348: 177-199. DOI:10.1017/S0022112097006551 [7] AUBRY N, GUYONNET R, LIMA R. Spatio-temporal analysis of complex signals:Theory and applications[J]. Journal of Statistical Physics, 1981, 64(3-4): 683-739. [8] HOLMES P, LUMLEY J L, BERKOOZ G. Turbulence, coherent structures, dynamical systems and symmetry[M]. Cambridge: Cambridge University Press, 1996: 68-100. [9] VOLKWEIN S.Proper orthogonal decomposition for nonlinear dynamical systems[EB/OL].Graz: University of Graz, 2005[2018-08-14] http://www.math.unikonstanz.de/numerik/personen/volkwein/PhDSchools/Volkwein_Part1.pdf. [10] NAKAKITA K, HABASHI W G, NADARAJAH S. Toward real-time aero-icing simulation using reduced order models[J]. Journal of Aircraft, 2010, 47(1): 96-115. DOI:10.2514/1.44077 [11] 葛宜元. 试验设计方法与Design-Expert软件应用[M]. 哈尔滨: 哈尔滨工业大学出版社, 2015: 2-3. GE Y Y. Experimental design method and Design-Expert software application[M]. Harbin: Harbin Institute of Technology Press, 2015: 2-3. (in Chinese) [12] MCKAY M D, BECKMAN R J, CONOVER W J. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technometrics, 1979, 21(2): 239-245. [13] MORRIS M D, MITCHELL T J. Exploratory designs for computational experiments[J]. Journal of Statistical Planning and Inference, 1995, 43(3): 381-402. DOI:10.1016/0378-3758(94)00035-T [14] 丁鹏, 陶文诠. 建立低阶模型的POD方法[J]. 工程热物理学报, 2009, 30(6): 1019-1021. DING P, TAO W Q. Reduced order modeling with the proper orthogonal decomposition[J]. Journal of Engineering Thermo Physics, 2009, 30(6): 1019-1021. DOI:10.3321/j.issn:0253-231X.2009.06.032 (in Chinese) [15] KERSCHEN G, GOLINVAL J, VAKAKIS A F, et al. The method of proper orthogonal de-composition for dynamical characterization and order reduction of mechanical sys-tems:An overview[J]. Nonlinear Dynamics, 2005, 41(1-3): 147-169. DOI:10.1007/s11071-005-2803-2 [16] SIROVICH L. Turbulence and the dynamics of coherent structures.Ⅰ-Coherent structures.Ⅱ-Symmetries and transformations.Ⅲ-Dynamics and scaling[J]. Quarterly of Applied Mathematics, 1987, 45(3): 561-571. DOI:10.1090/qam/1987-45-03 [17] KRIGE D G. A statistical approach to some basic mine valuation problems on the Witwatersrand[J]. Journal of the Southern African Institute of Mining and Metallurgy, 1951, 52(6): 119-139. [18] 韩忠华. Kriging模型及代理优化算法研究进展[J]. 航空学报, 2016, 37(11): 3197-3225. HAN Z H. Kriging surrogate model and its application to design optimization:A review of recent progress[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3197-3225. (in Chinese) [19] HÉLOÏSE B, FRANÇOIS M, WAGDI G H. FENSAP-ICE's three-dimensional in-flight ice accretion module:ICE3D[J]. Journal of Aircraft, 2003, 40(2): 239-247. DOI:10.2514/2.3113 [20] JEONG S, OBAYASHI S, YAMAMOTO K. Aerodynamic optimization design with kriging model[J]. Transactions of the Japan Society for Aeronautical and Space Sciences, 2005, 48(161): 161-168. DOI:10.2322/tjsass.48.161 [21] JECK R K.Icing Design Envelopes(14 CFR Parts 25 and 29, Appenddix C)Converted to a Distance-Based Format: DOT/FAA/AR-00/30[R].Washington, D.C.: FAA, 2002. [22] JOSEPH V R, HUNG Y, SUDJANTO A. Blind Kriging:A new method for developing metamodels[J]. Journal of Mechanical Design, 2008, 130(3): 350-353. [23] COUCKUYT I, FORRESTER A, GORISSEN D, et al. Blind Kriging:Implementation and performance analysis[J]. Advances in Engineering Software, 2012, 49: 1-13. DOI:10.1016/j.advengsoft.2012.03.002

#### 文章信息

LIU Teng, LI Dong, HUANG Ranran, ZHANG Zhenhui

Ice shape prediction method of aero-icing based on reduced order model

Journal of Beijing University of Aeronautics and Astronsutics, 2019, 45(5): 1033-1041
http://dx.doi.org/10.13700/j.bh.1001-5965.2018.0474