ADAPTIVE MULTI-FIDELITY POLYNOMIAL CHAOS-KRIGING MODEL-BASED EFFICIENT AERODYNAMIC DESIGN OPTIMIZATION METHOD
-
摘要: 多可信度代理模型已经成为提高基于代理模型的优化算法效率和可信度水平最有效的手段之一. 然而目前流行的co-Kriging和分层Kriging (HK)等多可信度代理模型泛化能力不足, 缺乏对高阶/高非线性建模问题的适应性, 难以广泛应用. 文章基于发展的自适应多可信度多项式混沌-Kriging (MF-PCK)代理模型, 在提高建模效率和对高阶/高非线性问题近似准确率的同时, 建立了基于该自适应MF-PCK模型的高效全局气动优化方法. 在发展的方法中, 提出了基于MF-PCK模型的新型变可信度期望改进加点方法, 使代理优化算法效率进一步提高. 为了验证发展方法的全面表现, 将其应用在经典的数值函数算例以及多个跨音速气动外形的确定性优化和稳健优化设计中, 并与基于Kriging和HK模型的代理优化算法进行了全面比较. 结果表明, 发展的新型多可信度全局气动优化方法其优化效率相对于基于Kriging和HK模型的优化效率显著提高, 结果更好也更加可靠, 并且稳健优化设计效率和结果也更符合工程应用需求, 证明了其相对于基于Kriging和HK模型的代理优化算法的显著优势.
-
关键词:
- 多可信度代理模型 /
- 全局优化 /
- 跨音速气动优化 /
- 多可信度多项式混沌-克里金 /
- 稳健设计
Abstract: The multi-fidelity surrogate-assisted global aerodynamic optimization method has become one of the most efficient means to improve the efficiency and reliability of surrogate-based optimization algorithms. However, currently popular multi-fidelity surrogate models, e.g., co-Kriging and hierarchical Kriging (HK) models, lack the generalization ability and the adaptability to complex engineering applications, so that it is very difficult to promote them to be widely utilized. In this paper, based on the developed adaptive multi-fidelity polynomial chaos-Kriging (MF-PCK) surrogate model, an efficient global aerodynamic optimization framework based on the novel adaptive MF-PCK model is proposed while improving the modeling efficiency and the global approximation accuracy for highly-nonlinear and high-order response problems. In the development framework, we proposed a new adaptive multi-fidelity polynomial chaos-Kriging-based variable-fidelity expected improvement (MF-PCK-VF-EI) infilling strategy that can adaptively select new samples of both low- and high-fidelity, which further improves the efficiency and ability of this optimization framework. The developed methodology is verified by one analytical function example and demonstrated by three engineering application problems, including the deterministic optimization and the robust design optimization of several transonic aerodynamic shapes, compared with the popular Kriging-based and HK-based optimization algorithms. The results show that the global aerodynamic optimization efficiency of the developed method based on the novel adaptive MF-PCK model is higher than that based on the HK model and more than twice that based on the Kriging model, meantime with better and more reliable results. More importantly, the efficiency and performance of robust aerodynamic design optimization (RADO) for transonic aerodynamic shapes are prone to satisfy the requirements for complex engineering applications, which proves a significant advantage of the developed method over the popular Kriging-based and HK-based optimization algorithms. -
引 言
随着人们对设计外形的性能表现和可信度要求的不断提高, 高可信度计算流体力学(computational fluid dynamics, CFD)数值模拟希望被应用到气动设计中, 这使得相同数量样本的CFD分析计算花费急剧增加, 因此急需发展更有效的代理优化方法[1-4]. 为此, 最流行的手段之一就是利用多可信度代理模型代替原来的单可信度代理模型以减少代理建模和优化过程所使用的高可信度样本数量. 多可信度代理模型通过使用大量且廉价的低可信度样本建模来反映高可信度模型趋势, 并利用少量的高可信度样本对该模型进行修正, 从而显著减小构造高可信度代理模型所需的计算花费(更少的高可信度样本数), 提高了代理建模和代理优化效率[5].
多可信度代理模型主要有3种构建方式[6], 即加法或乘法标度、混合标度以及空间映射. 加法标度相对于乘法标度能够更好地全局近似高可信度分析函数, 鲁棒性也更好. 典型的如多可信度多项式混沌展开(polynomial chaos expansion, PCE)[7]使用加法标度对低可信度稀疏PCE中重要项的恢复系数进行修正, 以减少所需的高可信度样本数. 其建模效率相对于单可信度稀疏PCE提高一倍多. 但多可信度PCE方法[7]缺乏合适的标准确定多项式矫正项项数, 对复杂的问题需要多次试探, 从而引入了额外的计算花费. 然而, 单一的加法或乘法标度缺乏对更多问题的适应性[6], 为此研究者又提出了混合标度法, 即同时使用加法和乘法标度修正低可信度模型. 在各种混合标度法中, 两种最流行的形式包括混合桥函数法[8]和类似co-Kriging[9]和分层Kriging (hierarchical Kriging, HK)[10]的贝叶斯估计方法[6]. 混合桥函数法指使用系数分别对乘法修正因子和加法修正因子进行加权.
相比于混合桥函数法, 基于贝叶斯理论的co-Kriging[9]和HK对许多问题表现出了良好的灵活性和鲁棒性, 被广泛用于航空航天系统优化[1]和稳健设计[11]中. co-Kriging和HK使用缩放因子和静态高斯随机过程分别作为乘法标度项和加法标度项对低可信度Kriging模型进行修正, 从而获得对高可信度函数的近似. 随后, 研究者分别将co-Kriging和HK推广到任意多层可信度数据建模和优化设计中, 优势更加明显. 然而对于复杂气动外形精细化设计问题, 仍然难以避免相关矩阵规模和条件数过大导致的训练花费过高和精度急剧下降问题. 更严重的是, Toal[12]和Zhao等[13]均发现: 当高、低可信度模型之间偏差较复杂(高非线性/高阶)或高可信度模型是高非线性/高阶时, 类似co-Kriging和HK等多可信度代理模型使用缩放因子(对低可信度Kriging)和高斯随机过程项难以准确近似这部分复杂(高非线性)偏差. 即使使用更多高可信度样本, 其建模精度也难以和直接使用高阶多项式作为趋势函数的多项式混沌-Kriging (polynomial chaos-Kriging, PC-Kriging)等相比, 缺乏对飞行器复杂气动设计等许多工程问题所需的建模效率和泛化能力.
不仅如此, 低可信度模型的选择对多可信度代理模型的建模效率和精度有显著影响, 越好的低可信度模型应该是其与高可信度模型间的偏差越简单或者相对于高可信度模型更低阶[14]. 然而, 越好的低可信度模型可能要求更高的计算花费. Peherstorfer等[6]在多可信度代理模型综述中指出低可信度模型可以通过3种方式获得, 即简化的低可信度模型、基于投影的低可信度模型以及数据拟合的低可信度模型. 简化的低可信度模型指通过专业知识对高可信度模型进行适度简化, 如在流场特性分析时, 可以采用从直接数值模拟(direct numerical simulation, DNS)到大涡模拟(large eddy simulation, LES)、雷诺平均(Reynolds average Navier-Stokes, RANS)、Euler方程再到势流方程求解等, 或者对于确定的物理方程如RANS, 可使用从密网格再到不同密度的稀网格等, 他们依次降低可信度和计算花费. 而基于投影的低可信度模型指通过数学的方法探测问题内在的结构, 比如通过在低维子空间投影得到的降阶模型等. 数据拟合类低可信度模型指构建输入到输出的低可信度映射或数学近似模型等, 如低阶回归或插值等. 不同的多可信度代理模型对低可信度模型的要求和适应性不同, 也显著影响了他们的建模效率和泛化能力, 越好的多可信度代理模型应该能适应广泛的低可信度模型并提供良好的预测能力.
多可信度代理模型泛化能力的不足也严重影响了代理优化算法的表现[4]. 为了改进基于多可信度代理模型的优化算法的表现, 目前主要包括两方面有潜力工作. 其一, 改进多可信度代理模型的泛化能力[12]. 如前所述, 当前流行的多可信度代理模型如co-Kriging, HK以及多可信度多项式混沌(multi-fidelity polynomial chaos expansion, MF-PCE)等模型在不同问题中表现差异明显. 在实际问题中, 如果使用不当他们建模和优化效果可能很差. 基于此, 文献[13,15]提出了一种自适应多可信度PC-Kriging (multi-fidelity polynomial chaos-Kriging, MF-PCK)代理模型对这一问题进行了有效解决, 并通过在复杂数值函数和气动数据建模问题中验证表明: 其相对于流行的通用Kriging, MF-PCE和HK等对高非线性/高阶响应的泛化能力更强和建模效率更高. 其二, 发展更有效的多可信度优化加点准则[16]. 目前针对多可信度代理优化的加点策略研究仍比较少. Zhang等[17]针对分层Kriging (HK)模型提出了变可信度期望改进(variable-fidelity expected improvement, VF-EI)加点策略. Huang等[18]提出通过最大化增广EI (augmented expected improvement, AEI)函数来选取高或低可信度样本点. Jiang等[19]发展的基于变可信度Kriging的变可信度置信下届(VF-LCB)准则等. 然而, 此类多可信度优化加点策略对更多工程实际问题的适应性以及效率并没有得到广泛验证, 尤其难以适用于一般多可信度代理模型如多可信度多项式混沌以及多可信度径向基函数等.
如上所述, 如何发展新的多可信度代理模型以改进他们的泛化能力, 以及改进基于多可信度代理模型的全局优化能力和拓宽应用范围是当前亟需解决的难题. 针对第一个难题, 本文首先引入了作者提出的自适应多可信度多项式混沌-Kriging (MF-PCK)代理模型[13]. 该MF-PCK模型通过结合PCE的全局近似特性以及Kriging的局部插值特点, 自适应选择最优的多项式基函数和构建最优的多可信度代理模型. 该方法已经通过多个复杂数值函数算例和跨音速气动数据建模应用进行了广泛验证, 结果均表明, 其相对于分层Kriging, PC-Kriging以及多可信度PCE等流行代理模型全局建模效率和泛化能力均显著提高, 尤其对高阶/高非线性等复杂响应的近似准确率提升显著, 大大拓展了多可信度代理模型的应用范围和适应性. 针对第二个难题, 本文发展了基于MF-PCK代理模型的高维全局气动优化设计新方法, 和基于MF-PCK代理模型的全局优化加点方法, 显著改进了多可信度代理优化算法的优化效率和优化能力. 为了验证新方法的表现, 本文将其应用在经典的数值算例, 以及跨音速气动外形确定性优化和复杂的气动稳健优化设计应用中. 结果均表明, 基于MF-PCK代理模型的优化方法相对于基于Kriging的代理优化算法效率提高超过一倍, 并且结果更好也更可靠. 该方法的成功验证为未来发展基于MF-PCK代理模型的大规模变量全局优化设计新方法奠定了基础[20], 也为解决当前多可信度代理优化算法适应范围窄的难题提供了新思路.
1. 自适应多可信度PC-Kriging代理模型
在气动特性数值分析中, CFD模拟根据所使用的网格数量, 控制方程及近似假设等不同可以提供多个可信度系列样本数据. 低可信度模型能提供相对于高可信度模型可参考的趋势, 并且越低可信度的模型其计算花费也通常显著更低. 利用数量充足但计算花费非常低的低可信度样本来训练趋势模型将大大减少多可信度代理模型的建模花费. 本文结合多可信度建模技术发展了一种更加高效的代理建模途径.
1.1 多可信度多项式混沌-Kriging模型近似
多可信度建模技术通过有技巧地利用充足但廉价的低可信度数据去辅助多可信度代理建模, 在提高代理建模精度的同时, 显著降低计算花费[6-7, 14, 21]. 为此, 提出一种全新的混合标度形式来矫正低可信度PC-Kriging模型[13], 即
$$ \begin{split} & {y_h}({\boldsymbol{X}}) = {\alpha _0}\underbrace {\left( {\sum\limits_{i \in A} {{\beta _i}{\psi _i}({\boldsymbol{X}})} + {z_l}({\boldsymbol{X}})} \right)}_{{{\hat y}_l}({\boldsymbol{X}})}+ \\ &\qquad \underbrace {\overbrace {\sum\limits_{{c_i} \in {A_c}} {{\alpha _{{c_i}}}{\psi _{{c_i}}}({\boldsymbol{X}})} }^{{y_c}({\boldsymbol{X}})} + {z_h}({\boldsymbol{X}})}_{C({\boldsymbol{X}})} \end{split} $$ (1) 其中
$ {\alpha }_{0} $ 为多乘标度因子, 其暗示了高、低可信度模型之间的相关性,$C\left(\boldsymbol{X}\right)$ 为增加标度项.${\widehat{y}}_{l}\left(\boldsymbol{X}\right)$ 代表低可信度PC-Kriging预测器. MF-PCK使用基于LAR的低可信度PC-Kriging作为预测模型, 即$$ {\hat y_l}({\boldsymbol{X}}) = \sum\limits_{j = 1}^M {{\beta _{{i_j}}}{\psi _{{i_j}}}({\boldsymbol{X}})} + {\boldsymbol{r}}_l^{\text{T}}{\boldsymbol{R}}_l^{ - 1}({{\boldsymbol{Y}}_l} - {\boldsymbol{\varPsi \beta }}) $$ (2) 式中
$ {\boldsymbol{r}}_{l} = {\left[R\right({\boldsymbol{X}}^{\left(1\right)},\boldsymbol{X}),R({\boldsymbol{X}}^{\left(2\right)},\boldsymbol{X}),\cdots ,R({\boldsymbol{X}}^{\left({N}_{l}\right)},\boldsymbol{X}\left)\right]}^{\mathrm{T}} $ 表示低可信度观测样本$ {\boldsymbol{S}}_{l} $ 与未知点$ \boldsymbol{X} $ 之间的相关向量, 他们的计算均需要代入低可信度模型相关的超参数$ {\boldsymbol{\theta }}_{l} $ 和$ {q}_{l} $ , 具体可参考文献[13, 22].在式(1)中,
$ {y}_{c}\left(\boldsymbol{X}\right) $ 是对低可信度PCE的矫正修正项. 多项式校正展开项表达如下$$ {y_c}({\boldsymbol{X}}) = \sum\limits_{{c_i} \in {A_c}} {{\alpha _{{c_i}}}{\psi _{{c_i}}}({\boldsymbol{X}})} $$ (3) 其中
$ {\alpha }_{{c}_{i}} $ 和$ {A}_{c} $ 分别表示校正展开多项式项系数和相应的多项式项索引集合. 式(3)中校正展开多项式项必须是低可信度PC-Kriging模型中多项式项的子集. 本文设计了一种新的基于嵌套留一交叉验证(leave-one-out cross-validation, LOOCV)的LAR自适应基函数选择方法(即Nested LOOCV-LAR), 详情可见1.2节.通过标准“universal Kriging”拟合过程即可得到MF-PCK在未知点的预测模型, 可参考前期工作[13], 得到MF-PCK预测模型在未知点
$ \boldsymbol{X} $ 的预测值为$$ \begin{split} & {{\hat y}_h}({\boldsymbol{X}}) = {\alpha _0}{{\hat y}_l}({\boldsymbol{X}}) + \sum\limits_{i = 1}^{{M_c}} {{\alpha _{{c_i}}}{\psi _{{c_i}}}({\boldsymbol{X}})} + \\ &\qquad {\boldsymbol{r}}_h^{\text{T}}({\boldsymbol{X}}){\boldsymbol{R}}_h^{ - 1}\left( {{{\boldsymbol{Y}}_h} - {\boldsymbol{F\alpha }}} \right) \end{split}$$ (4) 其中,
$ {\widehat{y}}_{l}\left(\boldsymbol{X}\right) $ 是低可信度PC-Kriging预测器, 如式(2)所示.$ {M}_{c} $ (${M}_{c}\leqslant M$ )表示校正展开项中多项式项数量. 缩放因子以及校正展开多项式系数向量$ \boldsymbol{\alpha } $ 可通过广义最小二乘方法估计得到, 即$$ {\boldsymbol{\alpha }} = {\left( {{{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{R}}_h^{ - 1}{\boldsymbol{F}}} \right)^{{{ - }}1}}{{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{R}}_h^{ - 1}{{\boldsymbol{Y}}_h} $$ (5) 相关函数
$ {\boldsymbol{R}}_{h} $ 中未知超参数$ {\boldsymbol{\theta }}_{h} $ 和$ {q}_{h} $ 通过最大化似然函数得到$$ {\left\{{\theta }_{h},{q}_{h}\right\}}_{ML} = \underset{{\theta }_{h},{q}_{h}}{\mathrm{arg}\mathrm{min}}\left({\sigma }_{{z}_{h}}^{2}{\left|{{\boldsymbol{R}}}_{h}\right|}^{1/{N}_{h}}\right)\text{, }{\theta }_{{h}_{i}} > 0,{q}_{h} > 0 $$ (6) 即最大似然估计(MLE)优化所得. 静态过程方差
$ {\sigma }_{{z}_{h}}^{2} $ 由$ {\widehat{\sigma }}_{{z}_{h}}^{2} = {\left({\boldsymbol{Y}}_{h}-\boldsymbol{F}\boldsymbol{\alpha }\right)}^{\mathrm{T}}{\boldsymbol{R}}_{h}^{-1}\left({\boldsymbol{Y}}_{h}-\boldsymbol{F}\boldsymbol{\alpha }\right)/{N}_{h} $ 估计. 然后, 通过代入估计的超参数${{{\boldsymbol{\theta}} }}_{h}, {{q}}_{{h}}$ 以及高可信度观测样本${{{\boldsymbol{S}}}}_{h}$ , 可以计算出相关矩阵${{{\boldsymbol{R}}}}_{h}$ 和相关向量${{{\boldsymbol{r}}}}_{h}$ . 进一步, 可以方便地得到MF-PCK预测模型在测试点${{\boldsymbol{X}}}$ 的估计均方误差(mean-squared error, MSE), 即$$ s_{{y_h}}^2\left( {\boldsymbol{X}} \right)\; = \sigma _{{z_h}}^2\left[ {1.0 - {\boldsymbol{r}}_h^{\text{T}}{\boldsymbol{R}}_h^{ - 1}{{\boldsymbol{r}}_h} + {\boldsymbol{\varGamma }}{{\left( {{{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{R}}_h^{ - 1}{\boldsymbol{F}}} \right)}^{{\text{ - }}1}}{{\boldsymbol{\varGamma }}^{\text{T}}}} \right] $$ (7) 其中
$\boldsymbol{\varGamma } = {\boldsymbol{r}}_{h}^{\mathrm{T}}{\boldsymbol{R}}_{h}^{-1}\boldsymbol{F}-{\boldsymbol{\varLambda }}^{\mathrm{T}}$ . MF-PCK方法计算相关矩阵$ {\boldsymbol{R}}_{h} $ 和$ {\boldsymbol{R}}_{l} $ 更简单也更直接. MF-PCK预测的准确率与相关函数$R({{{\boldsymbol{X}}}}^{\left(i\right)},{{{\boldsymbol{X}}}}^{\left(j\right)};{{{\boldsymbol{\theta}} }}_{l},{q}_{l})$ 和$R({{{\boldsymbol{X}}}}^{\left(i\right)},{{{\boldsymbol{X}}}}^{\left(j\right)};{{{\boldsymbol{\theta}} }}_{h},{q}_{h})$ 中超参数${{{\boldsymbol{\theta}} }}_{l},{q}_{l},{{{\boldsymbol{\theta}} }}_{h},{q}_{h}$ 的最优估计密切相关. 一些研究已经指出, 在超参数接近原点的位置, 最大似然函数通常是多峰的. 因此, 通过有效的超参数优化方法以避免其陷入局部最优是非常有必要的. 本文使用全局粒子群算法与基于梯度的优化器SNOPT相结合, 以得到全局最优的超参数$ {\boldsymbol{\theta }}_{l},{q}_{l},{\boldsymbol{\theta }}_{h},{q}_{h} $ . 优化过程中优化范围定义为$\left[{10}^{-4}{\boldsymbol{\theta }}_{{\rm{ini}}},{10}^{4}{\boldsymbol{\theta }}_{{\rm{ini}}}\right]$ ,${\boldsymbol{\theta }}_{{\rm{ini}}}$ 是超参数初始值.1.2 基于LOOCV-LAR的自适应基函数选择方法
为了构建最优的MF-PCK预测模型, 本节设计了一个基于嵌套留一交叉验证−最小角回归(nested LOOCV-based LAR, LOOCV-LAR)方法来选择最优的多项式项集合, 以建立最优的MF-PCK预测模型. 该方法详细步骤如下.
(1) 选择合适的试验设计方法, 并分别抽取
$ {N}_{h} $ 个高可信度训练样本${\boldsymbol{S}}_{h} = {\left\{{{{\boldsymbol{X}}}}^{\left(i\right)}\right\}}_{i = 1}^{{N}_{h}}$ 以及$ {N}_{l} $ 个低可信度训练样本$ {\boldsymbol{S}}_{l} = {\left\{{\boldsymbol{X}}^{\left(i\right)}\right\}}_{i = 1}^{{N}_{l}} $ . 分别使用高可信度模型和低可信度模型计算他们的响应, 即${\boldsymbol{Y}}_{h} = \left({y}_{h}\right({{X}}^{\left(1\right)}),{y}_{h}({{X}}^{\left(2\right)}),\cdots , {y}_{h}({{X}}^{\left({N}_{h}\right)}\left)\right)^{\mathrm{T}}$ 和${\boldsymbol{Y}}_{l} = {\left({y}_{l}\right({{X}}^{\left(1\right)}),{y}_{l}({{X}}^{\left(2\right)}),\cdots ,{y}_{l}({{X}}^{\left({N}_{l}\right)}\left)\right)}^{\mathrm{T}}$ . 通过使用经典截断格式在$ p $ 多项式阶截断获得候选多变量正交多项式基函数集合$ \left\{{\psi }_{1}\right(\boldsymbol{X}),{\psi }_{2}(\boldsymbol{X}),\cdots ,{\psi }_{{M}_{p}}(\boldsymbol{X}\left)\right\} $ .(2) 利用低可信度样本集合
$ {\boldsymbol{S}}_{l} $ 和它们的响应集合$ {\boldsymbol{Y}}_{l} $ , 运行LAR算法获得候选的稀疏多项式基函数集合. LAR算法选择基函数的相对顺序决定了基函数对输出向量的贡献大小[23], 即$\{{\psi }_{{r}_{1}}\left(\boldsymbol{X}\right),{\psi }_{{r}_{2}}\left(\boldsymbol{X}\right),\cdots , {\psi }_{{r}_{{M}_{L}}}\left(\boldsymbol{X}\right)\}$ , 其中${M}_{L} = {{\rm{min}}}({M}_{p},{N}_{l}-1)$ .(3) 利用低可信度样本集合
$ {\boldsymbol{S}}_{l} $ 以及相应的响应集合$ {\boldsymbol{Y}}_{l} $ , 从上一步排序的稀疏多项式项集合中, 使用前$ M $ ($ M = 1, 2, \cdots ,{M}_{L}) $ 个多项式项$\left\{{\psi }_{{r}_{1}}\left(\boldsymbol{X}\right),{\psi }_{{r}_{2}}\left(\boldsymbol{X}\right),\cdots , {\psi }_{{r}_{M}}\left(\boldsymbol{X}\right)\right\}$ 作为趋势函数$ \mathrm{建}\mathrm{立}{M}_{L} $ 个PC-Kriging代理模型$ \left\{{\widehat{y}}_{{l}_{M}}\right(\boldsymbol{X}\left)\right\} $ . 代理模型的基函数个数依次增加.(4) 对于
$ i = 1,\mathrm{ }2,\mathrm{ }\cdots ,{N}_{h} $ :①对每个低可信度PC-Kriging代理模型
$ {\widehat{y}}_{{l}_{M}}\left(\boldsymbol{X}\right) $ 使用高可信度样本集合${\boldsymbol{S}}_{h}^{(-i)} = {\boldsymbol{S}}_{h}/ \left\{{\boldsymbol{X}}^{\left(i\right)}\right\}$ 和相应的响应集合${\boldsymbol{Y}}_{h}^{(-i)} = {\boldsymbol{Y}}_{h}/ \left\{{y}_{h}\left({\boldsymbol{X}}^{\left(i\right)}\right)\right\}\mathrm{建}\mathrm{立}M$ 个MF-PCK模型$\left\{{\widehat{y}}_{{h}_{M,{M}_{c}}}^{(-i)}\left(\boldsymbol{X}\right)\right\}$ . 其中, 在建立每个MF-PCK模型的过程中, 高可信度校正展开项$ {y}_{c}^{\left({M}_{c}\right)}\left(\boldsymbol{X}\right) $ 通过(3)式选择得到, 其项数${M}_{c} = 1,2,\cdots ,{{\rm{min}}}\left(M,{N}_{h}\right)$ , 即在校正展开中相应的多项式基集合为$\{{\psi }_{{i}_{1}}\left(\boldsymbol{X}\right),{\psi }_{{i}_{2}}\left(\boldsymbol{X}\right),\cdots , {\psi }_{{i}_{{M}_{c}}}\left(\boldsymbol{X}\right)\} \left({M}_{c} = 1, 2, \cdots ,{{\rm{min}}}\left(M,{N}_{h}\right)\right)$ .②计算每个多可信度代理模型
$ {\widehat{y}}_{{h}_{M,{M}_{c}}}^{(-i)}\left(\boldsymbol{X}\right) $ 在每个点$ \left\{{\boldsymbol{X}}^{\left(i\right)}\right\} $ 处的预测残差${\delta }_{M,{M}_{c}}^{\left(i\right)} = {\widehat{y}}_{{h}_{M,{M}_{c}}}^{(-i)}\left({{{\boldsymbol{X}}}}^{\left(i\right)}\right)-{{y}}_{h}\left({{{\boldsymbol{X}}}}^{\left(i\right)}\right)$ .(5) 计算每个MF-PCK代理模型
$ {\widehat{y}}_{{h}_{M,{M}_{c}}}^{(-i)}\left(\boldsymbol{X}\right) $ 的残差总和${\varepsilon }_{M,{M}_{c}} = {{Err}_{LOO}}_{M,{M}_{c}} = \displaystyle\sum\limits _{i = 1}^{{N}_{h}}{\left({\delta }_{M,{M}_{c}}^{\left(i\right)}\right)}^{2}/{N}_{h}$ . 找到最小的残差${\varepsilon }^{\mathrm{*}} = {\varepsilon }_{{M}^{\mathrm{*}},{M}_{c}^{\mathrm{*}}} = {{\rm{min}}}\left\{{\varepsilon }_{M,{M}_{c}}\right\}(M = 1,\mathrm{ }2,\cdots ,{M}_{L}, {M}_{c} = 1,\mathrm{ }2,\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\cdots , {{\rm{min}}}\left({M},{N}_{h}\right))$ 以及相应的最优$ {M}^{\mathrm{*}} $ 和$ {M}_{c}^{\mathrm{*}} $ , 则可以得到最相关的多项式基集合$\{{\psi }_{{r}_{1}}\left(\boldsymbol{X}\right),{\psi }_{{r}_{2}}\left(\boldsymbol{X}\right),\cdots , {\psi }_{{r}_{{M}^{*}}}\left(\boldsymbol{X}\right)\}$ 以及相应的最优矫正展开项集合$\{{\psi }_{{r}_{1}}\left(\boldsymbol{X}\right), {\psi }_{{r}_{2}}\left(\boldsymbol{X}\right),\cdots , {\psi }_{{r}_{{M}_{c}^{*}}}\left(\boldsymbol{X}\right)\}$ . 图1描述了最相关的低可信度多项式集合和最优的校正多项式选择过程.![]()
图 1 自适应选择最相关的低可信度多项式集合和最优的矫正多项式集合图示说明Figure 1. Sketch map in selecting the optimal cardinalities of the low-fidelity PCE and corresponding correction expansion for the optimal MF-PCK model2. 基于MF-PCK模型的自适应加点优化方法
2.1 基于MF-PCK代理模型的优化流程
如上所述, 相对于单可信度代理模型, 为了相同的近似准确率, 多可信度代理模型能显著减少计算花费. 因此, 应用多可信度代理模型到气动优化中能使高、低精度的样本数据有机结合, 以平衡分别基于高、低可信度数据优化的差异, 减少迭代周期和总花费. 本文以第1节提出的自适应MF-PCK代理模型为基础, 发展基于变可信度模型的加点准则. 基于MF-PCK代理模型的优化流程图如图2所示.
![]()
图 2 基于自适应MF-PCK的多可信度气动优化流程图Figure 2. The flow chart of adaptive MF-PCK-based multi-fidelity aerodynamic optimization algorithm2.2 变可信度期望改进加点方法
变可信度期望改进(VF-EI)加点策略是基于基本的EI加点方法和分层Kriging方法发展而来的[17], 本文在此基础上发展了基于新型MF-PCK的VF-EI加点策略. 由第1节发展的MF-PCK预测器在任意点
$ \boldsymbol{X} $ 的预测值为$$ \begin{split} & {{\hat y}_h}({\boldsymbol{X}}) = {\alpha _0}{{\hat y}_l}({\boldsymbol{X}}) + \sum\limits_{i = 1}^{{M_c}} {{\alpha _{{c_i}}}{\psi _{{c_i}}}({\boldsymbol{X}})} + \\ &\qquad {\boldsymbol{r}}_h^{\text{T}}({\boldsymbol{X}}){\boldsymbol{R}}_h^{ - 1}\left( {{{\boldsymbol{Y}}_h} - {\boldsymbol{F\alpha }}} \right) \end{split} $$ (8) 其中
${\widehat{{y}}}_{l}\left(\boldsymbol{X}\right)$ 为低可信度PC-Kriging预测器,$ \boldsymbol{F} $ 为高可信度样本的低可信度预测值集合$$ {\boldsymbol{F}} = {\left[ {{{\hat y}_l}({{{X}}^{(1)}}),{{\hat y}_l}({{{X}}^{(2)}}), \cdots ,{{\hat y}_l}({{{X}}^{({N_h})}})} \right]^{\text{T}}} $$ (9) 则在任意点的均方差预测(MSE)为
$$ \left.\begin{split} & s_h^2({\boldsymbol{X}}) = \sigma _{{z_h}}^2\left[ {1.0 - {\boldsymbol{r}}_h^{\text{T}}{\boldsymbol{R}}_h^{ - 1}{{\boldsymbol{r}}_h} + {\boldsymbol{\varPi}} {{\left( {{{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{R}}_h^{ - 1}{\boldsymbol{F}}} \right)}^{{{ - }}1}}{{\boldsymbol{\varPi}} ^{\text{T}}}} \right] \\ &{\boldsymbol{\varPi}} = {{\boldsymbol{r}}_h^{\text{T}}{\boldsymbol{R}}_h^{ - 1}{\boldsymbol{F}} - {{\boldsymbol{\varLambda }}^{\text{T}}}} \end{split}\right\} $$ (10) 在VF-EI中, MF-PCK预测的不确定性被定义为
$$ {s^2}({\boldsymbol{X}},l) = \left\{ \begin{gathered} \alpha _0^2s_l^2({\boldsymbol{X}}),\;\;\;\;\;\;\;l = 1\;{\text{for}}\;{\text{LF}}\;{\text{model}} \\ s_h^2({\boldsymbol{X}}),\;\;\;\;\;\;\;\;\;\;l = 2\;\;{\text{for}}\;{\text{HF}}\;{\text{model}} \\ \end{gathered} \right. $$ (11) 其中
$ {\alpha }_{0}^{2}{s}_{l}^{2} $ 为低可信度PC-Kriging预测的不确定性, 则相对于上式定义的双可信度预测器, VF-EI可定义为$$ {I_{vf}}({\boldsymbol{X}},l) = \max \left( {{y_{h,\min }} - {y_h}({\boldsymbol{X}}),0} \right) $$ (12) 相对于基于HK定义的VF-EI加点方法[17], 基于MF-PCK的VF-EI方法能提供更好的全局趋势预测能力和局部空间预测精度. 这是因为分层Kriging和co-Kriging等模型仅仅使用缩放因子和高斯随机过程去建模高、低可信度模型之间的差异部分, 而广泛的工程问题中高、低可信度模型之间的差异是非常复杂和高非线性的, 这使得他们难以提供准确的趋势预测模型, 也就导致了他们的均方差预测无法准确暗示误差最大位置. 而MF-PCK模型使用低可信度PC-Kriging和高可信度矫正展开多项式可以准确地近似高可信度函数趋势部分, 高斯随机过程可以精确模拟剩余残差, 使得MF-PCK均方差预测和预测误差一致性更好, 近似高可信度模型的准确率也显著更高[13], 更多详细的分析和验证可参考课题组相关的文献[13]. 通过相似的方式[24], 可以得到目标函数改进的期望为
$$ E({I_{vf}}) = \left\{ \begin{split} & \left( {{y_{h,\min }} - {f_h}({\boldsymbol{X}})} \right){\boldsymbol{\varTheta}} \left( {\frac{{{y_{h,\min }} - {f_h}({\boldsymbol{X}})}}{{s({\boldsymbol{X}},l)}}} \right) +\\ &\qquad s({\boldsymbol{X}},l)\varphi \left( {\frac{{{y_{h,\min }} - {f_h}({\boldsymbol{X}})}}{{s({\boldsymbol{X}},l)}}} \right),\;\;\;{\text{if}}\;\;s({\boldsymbol{X}}) > 0 \\ & 0,\qquad{\text{if}}\;\;s({\boldsymbol{X}}) = 0 \end{split} \right. $$ (13) 然后在VF-EI加点过程中的子优化就为寻找合适的加点位置以及可信度水平使VF-EI最大, 即
$$ \left.\begin{split} & {\rm{find}}\;\;{{\boldsymbol{X}}^*} \in [{{\boldsymbol{X}}_l},{{\boldsymbol{X}}_u}],l = 1,2 \\ & {\rm{max}}\;\;E({I_{vf}}({\boldsymbol{X}},l)) \\ & {\rm{s.t.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \;{\kern 1pt} {g_j}({\boldsymbol{X}},l) \geqslant 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \forall j = 1,2, \cdots ,{m_c} \end{split}\right\} $$ (14) 其中
$ {\mathrm{g}}_{j} $ 表示约束函数,$ {m}_{c} $ 为约束函数个数.3. 优化应用验证
为了检测提出的MF-PCK代理模型在优化设计中的表现, 本节分别使用一个典型函数优化算例和两个气动优化应用进行测试. 基于MF-PCK代理模型的优化过程使用VF-EI加点策略(MF-PCK_VF-EI), 并与使用基于Kriging的EI (Kriging_EI)加点算法和使用基于分层Kriging (HK)的VF-EI (HK_VF-EI)加点算法进行对比, 3种优化方法均使用全局加强学习粒子群(comprehensive learning particle swarm optimizer, CLPSO)算法.
3.1 函数优化测试
测试函数使用Branin函数, 其表达式为
$$ \begin{split} & {f_h}({\boldsymbol{X}}) = {\left({x_2} - \frac{{5.1}}{{4{\text{π} ^2}}}x_1^2 + \frac{5}{\text{π} }{x_1} - 6\right)^2} + \\ &\qquad 10\left(1 - \frac{1}{{8\text{π} }}\right)\cos {x_1} + 10,\; \\ &\qquad{x_1} \in [ - 5,10],\;{x_2} \in [0,15] \end{split} $$ (15) 相应的低可信度表达为
$$ {f_l}({\boldsymbol{X}}) = {f_h}({\boldsymbol{X}}) - 1.5{({x_2} - \frac{{5.1}}{{4{\text{π} ^2}}}x_1^2 + \frac{5}{\text{π} }{x_1} - 6)^2} $$ (16) Branin函数的函数值分布如图3所示, 图中红色实心点表示3个最小值位置, 最小值均为0.397887, 最小值位置分别为
$\left(-\mathrm{\text{π} },12.275\right),\mathrm{ }(\mathrm{\text{π} },2.275)$ 以及$(9.42\;478, \mathrm{ }2.475)$ . 对应Branin函数的低可信度函数如图4所示. 对基于MF-PCK_VF-EI和HK_VF-EI的多可信度优化过程, 通过拉丁超立方抽样(Latin hypercube sampling, LHS)分别选取5个高可信度样本点和50个低可信度样本点作为初始样本点, 而基于Kriging_EI的代理优化过程选取10个高可信度样本点作为初始样本. 3种方法均使用全局加强学习粒子群(CLPSO)优化算法. 当两种方法分别添加50个高可信度样本点时, 结果均接近收敛, 高可信度样本点对比分布如图5 ~ 图7所示. 基于MF-PCK_VF-EI和HK_VF-EI的方法在3个极值点附近均添加了更多样本点, 而基于Kriging_EI的方法添加样本点的位置仍然变化较大. 图8给出了函数值收敛过程对比, 明显可以看出, 基于MF-PCK_VF-EI的方法优化过程收敛最快, 其次是基于HK_VF-EI的方法, 并且两种多可信度代理优化方法均较基于Kriging_EI的优化算法收敛更快. 基于MF-PCK_VF-EI的优化方法在添加了20个高可信度样本以后就已经收敛, 而基于HK_VF-EI和Kriging_EI的优化方法则分别需要添加37和40个高可信度样本以后才接近收敛. 图9给出了每次加点位置处预测误差对比曲线. 可以看出, 基于MF-PCK_VF-EI方法的预测误差明显小于基于HK_VF-EI和Kriging_EI方法的预测误差, 并且较早收敛; 而基于HK_VF-EI和Kriging_EI的方法都需要添加到40个样本以后, 预测误差才不再明显波动. 表1给出了3种方法找到的最优位置和最优值对比, 可以看出, 在相同较少的计算花费下, 基于MF-PCK的方法找到的最优值低于基于HK_VF-EI和Kriging_EI的方法找到的最优值, 他们找到的位置也有较大区别.![]()
图 5 基于Kriging_EI加点过程Figure 5. The infilling process of Kriging_EI optimization method![]()
图 7 基于MF-PCK_VF-EI的加点过程Figure 7. The infilling process of MF-PCK_VF-EI optimization method![]()
图 9 代理模型预测误差收敛过程对比Figure 9. The convergence history of the surrogate prediction errors表 1 3种优化方法结果对比Table 1. Comparison of optimization results using three methodsOptimization method The optimal
locationsThe optimal values Prediction errors/10−4 Kriging_EI [−3.1702579,12.2190345] 0.4174474 8.953 HK_VF-EI [3.1912351, 2.3199426] 0.4166618 7.642 MF-PCK_VF-EI [9.4185951, 2.4763138] 0.3981135 3.692 3.2 跨音速气动优化应用
3.2.1 跨音速RAE2822翼型优化设计
首先测试了针对RAE2822翼型的不同网格水平的计算结果. 网格为C型拓扑, 远场距离为50倍翼型弦长, 湍流模型均使用两方程SST. 不同网格量网格近场分布如图10所示. 远场网格分布如图11所示. 选取的典型计算状态为:
$ Ma = 0.734 $ ,$Re = 6.5\times {10}^{6}$ , 以及${{C}}_{l} = 0.824$ ($\alpha = {2.79}^{\circ }$ ). 图12给出了数值模拟结果与试验值压力分布对比. 可以看出, 随着网格密度增加, 数值模拟结果与试验结果吻合越来越好, 对激波位置和强度的捕捉也更为准确. 表2给出了数值计算结果与试验气动力系数的对比. 结果显示, 从稀网格($ {L}_{1} $ 水平)增加到中等网格($ {L}_{3} $ 水平), 网格量增加了接近20.25倍, 阻力系数误差降低了17.22%. 继续增加网格量到密网格($ {L}_{5} $ 水平), 网格量增加了49.5倍, 阻力系数降低了18.80%, 而与$ {L}_{3} $ 水平仅仅相差了1.58%, 压力系数分布接近收敛. 因此本文分别使用$ {L}_{1} $ 水平和$ {L}_{3} $ 水平为低可信度和高可信度样本CFD分析网格, 使用$ {L}_{5} $ 水平网格作为优化外形气动特性检验网格, 以满足高效气动优化设计需求.![]()
图 10 针对RAE2822翼型的不同可信度水平计算网格Figure 10. Different levels of fidelity grids around RAE2822 airfoil![]()
图 12 不同网格量计算压力分布与试验压力分布对比Figure 12. Comparison of pressure distributions among computational and experimental results表 2 RAE2822网格收敛结果对比Table 2. Comparison of computational results of different levels of gridsGrid level Grid number $ {C}_{l} $ $ {C}_{d} $ $ {C}_{m} $ α/(°) $ {L}_{1} $ $4.0\times10^3$ 0.824 0.024297 −0.096164 3.022 $ {L}_{2} $ $ 2.1\times {10}^{4} $ 0.824 0.021261 −0.092284 2.973 $ {L}_{3} $ $ 8.1\times {10}^{4} $ 0.824 0.020112 −0.092868 2.917 $ {L}_{4} $ $1.27\times10^5$ 0.824 0.019960 −0.096062 2.869 $ {L}_{5} $ $1.98\times10^5$ 0.824 0.019729 −0.092778 2.806 为了对比说明, 本节基于单可信度Kriging模型的初始样本为
$ {N}_{h} = 40 $ , 即计算花费为$ 40 t $ . 而基于MF-PCK模型的方法使用课题组提出的基于留一交叉验证-Voronoi-最大尺度距离(LOOCV-Voronoi-MSD)的序列抽样方法[13]从$ {N}_{h} = 1 $ 和$ {N}_{l} = 10 $ 开始自适应选样直到获得相同计算花费($ 40 t $ )的高、低可信度样本分布和样本数量, 即使用$ {N}_{h} = 15 $ 个高可信度样本和$ {N}_{l} = 500 $ 个低可信度样本. 其中$ t $ 表示一次高可信度CFD计算分析时间, 高、低可信度数值分析均使用RANS数值求解器求解. 设计状态为$ Ma = 0.734,{C}_{l} = 0.824,Re = 6.5\times {10}^{6} $ . 设计约束包括升力系数、力矩系数和面积均不减小, 优化模型为$$ \left.\begin{split} & {\rm{find}}\;\;\;{\boldsymbol{X}} \in {{\boldsymbol{X}}^d} \\ & {\rm{min}}\;\;\;{C_d} \\ & {\rm{ s.t.}}\;\;\;\;\;{C_l} = 0.824,\;{C_m} \geqslant - 0.092, \\ &\qquad {{\rm{area}}_{{\rm{optimum}}}} \geqslant {{\rm{area}}_{{\rm{baseline}}}} \end{split}\right\} $$ (17) 本算例使用SST湍流模型进行气动分析, 初始翼型阻力为0.02171. 采用类函数/形函数(class/shape transformation, CST)参数化方法控制翼型变形, 上下表面各选取12个设计变量, 共24个变量, 初始设计范围如图13所示. 3种方法均使用全局加强学习粒子群(CLPSO)优化算法. 基于MF-PCK_VF-EI的优化流程如图2所示, 迭代加点70次. 图14给出了3种方法阻力系数收敛曲线对比, 可以看出, 相比于基于Kriging_EI和HK_VF-EI的两种优化方法, 基于MF-PCK_VF-EI的代理优化算法获得了最快的收敛速度以及最低的阻力适应值, 在迭代到第30次高可信度加点时, 阻力系数已经收敛; 而基于HK_VF-EI的优化算法则需要迭代加点到50个高可信度加点时才逐渐收敛; 基于Kriging_EI的代理优化方法收敛速度最慢, 直到第70次迭代加点时阻力系数仍在少量下降. 从图15可以找出原因, MF-PCK代理模型在迭代过程中持续拥有更准确的预测精度, 从第30次高可信度加点以后, 加点位置处阻力系数预测误差在0.0002以内, 而HK和Kriging代理模型阻力系数预测误差仍然在0.0006较大范围以内波动. MF-PCK模型使用了较多的低可信度样本, 可以更准确地捕捉整个设计空间的特点. 在Kriging_EI加点过程中, 由于EI搜索兼具了提高全局近似准确率的特性, 造成了加点位置波动较大. 当加点样本较少时, Kriging代理模型本身的近似不足, 进而造成最优位置以及近似误差的波动. 随着基于Kriging_EI的加点过程进行, 优化目标值和代理预测误差均逐渐收敛, 但优化目标的收敛却慢于代理预测误差的收敛, 这与EI加点过程后期收敛性变差密切相关, 当代理预测误差和优化目标均收敛时, 该优化过程趋于收敛, 得到收敛外形. 表3给出了基于3种代理模型的优化算法在接近收敛时的总计算时间和气动力系数对比. 基于MF-PCK_VF-EI的代理优化方法在第30次高可信度加点时阻力系数已经收敛, 收敛阻力系数值为111.82 counts, 总计算时间(初始样本和新添加样本的CFD分析时间)约等于
$ 70 t $ ; 而基于HK_VF-EI和Kriging_EI的代理的优化方法分别在第50次和第70次加点时阻力系数在接近收敛, 收敛的阻力系数值分别为112.82 counts和113.75 counts, 总计算时间分别约等于$ 90 t $ 和$ 110 t $ . 因此, 基于MF-PCK的代理优化方法相对于HK_VF-EI和Kriging_EI等两种经典代理优化方法在更短时间内获得了阻力值更低的翼型, 效率和可靠性均显著提升. 此外, 因为在基于MF-PCK_VF-EI的优化过程中, 子迭代自适应选择和分析低可信度样本点的计算花费非常低, 其相对于添加高可信度样本点的花费可忽略不计, 因此本文忽略了加点过程中低可信度样本的计算花费. 图16给出了优化翼型和RAE2822翼型外形对比, 可以看出基于MF-PCK_VF-EI的方法优化的翼型较基于Kriging_EI方法优化的翼型最大厚度位置略有后移, 弯度略有增加. 图17给出了4个翼型的压力分布对比, 可以看出, 基于MF-PCK_VF-EI的方法优化的翼型激波几乎完全抹平.![]()
图 15 阻力系数代理模型预测误差收敛过程Figure 15. The convergence history of surrogate prediction errors of drag coefficients表 3 RAE2822翼型优化设计结果对比Table 3. Comparison of optimization results of RAE2822 airfoilOptimization methods $ {C}_{l} $ ${C}_{d}/10^{-3}$ Areas ${C}_{m}/10^{-3}$ Total time of CFD evaluations Kriging_EI 0.824 11.375 0.0779 −91.98 $ 70{N}_{h} + 40{N}_{h} = 110 t $ HK_VF-EI 0.824 11.282 0.0779 −91.75 $ 50{N}_{h} + 15{N}_{h} + 500{N}_{l}\approx 90 t $ MF-PCK_VF-EI 0.824 11.182 0.0779 −91.71 $ 30{\mathit{N}}_{\mathit{h}} + 15{\mathit{N}}_{\mathit{h}} + 500{\mathit{N}}_{\mathit{l}}\approx 70\mathit{t} $ ![]()
图 17 优化翼型压力分布对比Figure 17. Comparison of pressure distributions of optimized and initial airfoils3.2.2 跨音速RAE2822翼型稳健优化设计
然而当考察单点优化得到的翼型在非设计点表现时, 结果如图18所示. 单点优化翼型在非设计点阻力剧烈增加, 其气动表现对马赫数、升力系数等状态参数波动非常敏感, 因而难以被工程使用. 针对这个问题, 稳健优化是最有希望的解决途径[25-27]. 然而, 稳健优化是一个典型的双循环过程[26, 28], 外循环为优化设计循环, 内循环在迭代的每一步对每个候选设计(粒子)进行不确定分析, 得出量化的不确定量, 如均值和方差等. 这使得即使在内循环使用高效的不确定分析方法, 如PCE, 但依然不可避免的要求远远更多的计算花费, 远超过基于代理的确定性优化[11, 29]. 针对这个问题, 本文构建关于设计变量与随机变量的联合代理模型, 使得在搜索过程中通过调用此代理模型对每个候选设计进行不确定分析, 从而大大减少了计算花费[30-31]. 定义上述问题的稳健优化模型为
![]()
图 18 单点优化设计翼型在非设计点气动特性对比Figure 18. Comparison of aerodynamic characteristics of optimized airfoil on off-design points$$ \left.\begin{split} & {\rm{condition}}:\;\;{{{Re}}} = 6.5 \times {10^6},\;{C_l} = 0.824, \\ & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;M{a_\infty } \in U(0.7,0.75)\;\;\;\;\;\;\;\;\;\;\;\;\; \\ & {\rm{find}}:\;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{X}} \in {R^n},\;{\boldsymbol{X}} \in [{{\boldsymbol{X}}_L},{{\boldsymbol{X}}_U}] \\ & {\rm{minimize}}:\;\;\;\;\;\left\{ \begin{gathered} {u_{_{{C_d}}}}(Ma) \\ {\sigma _{{C_d}}}(Ma) \\ \end{gathered} \right. \\ & {\rm{subject}}\;\;{\rm{to}}:\;\;\;Pr ({C_m} \leqslant {C_{m0}}) \leqslant {P_0} \\ & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{area}} \geqslant {{\rm{area}}_{{\rm{rae}}2822}} \end{split}\right\} $$ (18) 对该优化问题首先构建马赫数
$ Ma\in \left[\mathrm{0.7,0.75}\right] $ 与设计变量${x}_{i}\in \left[-\mathrm{0.03,0.03}\right](i = \mathrm{1,2},\cdots, 24)$ 共25个变量的联合代理模型$ S(Ma,{x}_{1},{x}_{2},\cdots ,{x}_{24}) $ , 其中设计变量设计范围如图13所示. 在联合代理模型$\widehat{y}(\boldsymbol{X},\boldsymbol{\varXi })$ 构建后, 在迭代的每一步, 每个候选设计的表现不确定量化可通过下式估算$$ \left.\begin{split} & {\mu _y}({\boldsymbol{X}}) = \frac{1}{N}\sum\limits_{i = 1}^N {\hat y({\boldsymbol{X}},{{\boldsymbol{\varXi }}_i})} \\ & {\sigma _f}({\boldsymbol{X}}) = \sqrt {\frac{1}{{N - 1}}\sum\limits_{i = 1}^N {{{(\hat y({\boldsymbol{X}},{{\boldsymbol{\varXi }}_i}) - {\mu _f}({\boldsymbol{X}}))}^2}} } \end{split}\right\} $$ (19) $$ \left.\begin{split} & {P_g} = \frac{1}{N}\sum\limits_{i = 1}^N {{I^{(i)}}} \;\;\; \\ & {I^{(i)}} = I({\boldsymbol{X}},{{\boldsymbol{\varXi }}^{(i)}}) = \\ &\qquad \left\{ \begin{gathered} 1,\;\;{\rm{if}}\;g({\boldsymbol{X}},{{\boldsymbol{\varXi }}^{(i)}}) \leqslant 0 \\ 0,\;\;{\rm{if}}\;g({\boldsymbol{X}},{{\boldsymbol{\varXi }}^{(i)}}) > 0\; \\ \end{gathered} \right. \;\;\;(i = 1,2, \cdots N) \end{split}\right\} $$ (20) 其中
$\boldsymbol{\varXi }$ 代表了不确定因素, 如Ma.通过LHS抽样选取100个混合变量高可信度样本以及1000个混合变量低可信度样本, 构建联合代理模型. 此时, MF-PCK预测力矩系数均值和标准差RRMSE分别为0.00364和0.01935, 完全可以满足优化设计使用. 然后使用全局加强学习粒子群(CLPSO)优化算法优化迭代1000次, 得到如图19所示稳健优化翼型. 可以发现, 相比于单点优化翼型, 最大厚度位置后移, 前缘半径变小. 图20给出了稳健设计翼型与单点优化翼型及初始翼型阻力发散特性对比. 可以看出, 稳健优化翼型阻力发散马赫数超过0.75, 高于其他两个翼型, 并且在马赫数0.70 ~ 0.75范围内阻力系数变化更加稳健. 图21给出了稳健设计翼型与其他两个翼型在不同马赫数下压力分布对比. 在0.734马赫数下, 稳健优化翼型上表面出现一个双弱激波, 在0.75马赫数下, 稳健优化翼型也出现了一个双弱激波, 强度均明显弱于单点设计翼型的强激波. 因此, 通过稳健设计优化, RAE2822翼型上表面出现了稳定的双弱激波, 从而提高了对马赫数的稳健性. 然而, 针对其他气动外形设计算例, 双弱激波是否具有对马赫数等的稳健性, 仍然值得更多详细研究.
![]()
图 20 优化翼型阻力发散特性对比Figure 20. Comparison of drag-divergence characteristics of optimized and initial airfoils![]()
图 21 优化翼型与RAE2822翼型在不同状态下压力分布对比Figure 21. Comparison of pressure distributions of optimized and initial airfoils under different Mach numbers3.2.3 跨音速ONERA M6机翼优化设计
跨音速ONERA M6机翼经常被用来测试数值计算程序或者优化算法的可靠性. 本算例考虑的计算状态为
$Ma = 0.839\;5,Re = 1.172\times {10}^{7},\alpha = {3.06}^{\circ }$ . 该状态下, ONERA M6机翼上表面有较强激波. 首先使用3套不同网格量计算网格验证计算程序的可靠性, 计算网格分布如图22所示. 图23分别给出了41%和63%位置处剖面压力分布对比, 可以看出, 640 W网格和试验压力分布吻合良好, 而320 W网格计算结果接近640 W网格计算结果, 30 W网格计算结果能捕捉到激波基本位置. 正如第1节指出, 低可信度计算网格的选择应该以能较为准确地捕捉到高可信度计算的趋势为准. 因此, 在本算例中, 高可信度气动分析使用网格量为320 W的计算网格, 而低可信度计算使用网格量为30 W的网格, 两者在相同计算资源下, 计算时间相差10.67倍. 本算例均使用SST湍流模型进行计算, 高可信度计算初始机翼阻力为113.97 counts, 升力系数为0.1670. 设计问题定义为在机翼各个剖面最大厚度和机翼升力系数不减小的情况下最小化机翼阻力系数, 其优化数学模型为![]()
图 23 不同计算网格量下计算压力分布与试验压力分布对比Figure 23. Comparison of computational and experimental pressure distributions under different number of grids$$ \left.\begin{split} & {\rm{condition}}:\;\;{{Ma}} = 0.839\;5,{{Re}} = 1.172 \times {10^7},\\ &\qquad\qquad\qquad{{\alpha}} = {3.06^\circ} \\ & {\rm{find}}:\;\;\;\;\;\;\;\;{\boldsymbol{X}} \in {R^n},\;{\boldsymbol{X}} \in [{{\boldsymbol{X}}_L},{{\boldsymbol{X}}_U}] \\ & {\rm{min}}:\;\;\;\;\;\;\;\;\;{C_D} \\ & {\rm{s.t.}}\;\;\;\;\;\;\;\;\;\;\;{C_l} \geqslant {C_{{l_0}}},\;{{\rm{thick}}_{{\rm{optimum}}}} \geqslant {{\rm{thick}}_{{\rm{baseline}}}} \end{split}\right\} $$ (21) 采用自由变形(free form deformation, FFD)方法对机翼进行参数化, 如图24所示, 选择4个剖面对机翼进行变形, 每个剖面定义了12个控制点, 通过控制每个控制点沿z方向移动从而产生新外形, 机翼剖面不扭转. 针对基于MF-PCK_VF-EI的代理优化方法, 采用LHS产生
${N}_{{\rm{ini}},h} = d\times 2 = 96$ 个初始高可信度样本, 以及${N}_{{\rm{ini}},l} = d\times 20 = 960$ 个初始低可信度样本,采用VF-EI加点方法进行优化. 两种方法均使用全局加强学习粒子群(CLPSO)优化算法. 作为对比, 针对基于Kriging代理模型的优化方法使用LHS产生${N}_{{\rm{ini}}} = d\times 5 = 240$ 个初始高可信度样本, 采用EI加点优化方法, 其中$ d $ 为设计变量数. 图25给出了两种方法阻力收敛历程对比. 结果显示, 基于MF-PCK的代理优化方法随着高可信度样本填充, 具有更快的收敛速度, 到78次高可信度加点时阻力已经收敛到98.07 counts, 阻力下降了15.9 counts; 而基于Kriging_EI的代理优化方法收敛很慢, 直到200次加点时阻力才收敛到99.04 counts. 因此, 基于提出的多可信度代理模型采用更少的计算花费获得了更好的收敛解, 大大提高了代理优化收敛速度以及获得更优解的效率(优化效率提高了一倍多). 图26给出了初始M6机翼和(基于MF-PCK)代理优化的机翼上表面压力云图对比. 可以看出, 优化的机翼基本消除了原始机翼上表面的两道强激波, 仅仅在局部位置存在非常弱的激波. 可以发现到接近收敛时, 每一次迭代阻力值下降很小. 这是因为在无激波外形附近, 阻力值变化接近平坦, 随机搜索算法效率很低. 图27给出了4个剖面位置处的压力分布和剖面翼型对比, 强激波已经消除. 设计结果再次验证了提出的方法相对于基于Kriging_EI的代理优化方法的有效性和巨大优势.![]()
图 25 两种方法阻力收敛历程对比Figure 25. The convergence history of different optimization algorithms![]()
图 26 初始机翼和优化后机翼(MF-PCK_VF-EI)上表面压力云图对比Figure 26. Comparison of pressure distribution contours of upper surface of M6 wing using the MF-PCK_VF-EI method![]()
图 27 初始机翼和优化后的机翼(MF-PCK_VF-EI)各剖面压力分布和翼型对比Figure 27. Comparison of pressure distributions and shapes of initial and optimized wing sections4. 总结
本文针对基于代理模型的全局气动优化方法效率低下和CFD可信度要求不断提高带来的优化设计难题, 提出了基于新型自适应多可信度多项式混沌-Kriging (MF-PCK)代理模型的高效全局气动优化方法和框架. 本文详细推导了建模过程及构建步骤, 并陈述了优化框架中的关键模块和技术难点. 通过使用相同标度形式再次矫正MF-PCK模型, 该方法也可用于3层及以上多水平多可信度气动数据建模和优化设计中. 同时本文提出了基于MF-PCK模型的变可信度期望改进加点优化方法, 其在每次加点迭代过程中能自适应地选择高可信度或低可信度样本添加到样本库中, 提高优化过程中气动分析效率, 从而进一步提升多可信度全局气动优化效率和能力. 最后本文将发展的框架应用到RAE2822翼型和M6机翼的单点确定性设计和RAE2822翼型的稳健优化设计中, 并与基于HK_VF-EI和Kriging_EI的代理优化方法进行了全面比较. 基于本文气动优化算例的结果显示, 基于新型MF-PCK代理模型的优化方法其优化效率相对于基于HK_VF-EI和Kriging_EI的优化效率显著提高, 结果更好也更加可靠, 并且稳健优化设计效率和结果也更加符合实际工程应用对气动优化方法的要求, 从而广泛验证了发展的新方法相对于当前基于HK_VF-EI和Kriging_EI的代理优化算法的显著优势和巨大潜力. 此外, 虽然MF-PCK模型相比HK模型等具有更强的泛化能力, 但训练时间也更长, 未来可通过发展基于多核并行的自适应LOOCV-LAR基函数选择方法来有效解决. 同时当前流行的多可信度代理模型如Co-Kriging, HK以及MF-PCK等均遭遇维数灾难难题, 未来需要发展更有效的高维空间代理建模技术, 如降维代理建模方法等. 最后基于MF-PCK代理模型的优化算法亦缺乏对广泛工程问题的验证, 未来需要结合具体的复杂应用来探究它的实际的表现和可能的改进.
-
![]()
图 1 自适应选择最相关的低可信度多项式集合和最优的矫正多项式集合图示说明
Figure 1. Sketch map in selecting the optimal cardinalities of the low-fidelity PCE and corresponding correction expansion for the optimal MF-PCK model
![]()
图 2 基于自适应MF-PCK的多可信度气动优化流程图
Figure 2. The flow chart of adaptive MF-PCK-based multi-fidelity aerodynamic optimization algorithm
![]()
图 5 基于Kriging_EI加点过程
Figure 5. The infilling process of Kriging_EI optimization method
![]()
图 7 基于MF-PCK_VF-EI的加点过程
Figure 7. The infilling process of MF-PCK_VF-EI optimization method
![]()
图 9 代理模型预测误差收敛过程对比
Figure 9. The convergence history of the surrogate prediction errors
![]()
图 10 针对RAE2822翼型的不同可信度水平计算网格
Figure 10. Different levels of fidelity grids around RAE2822 airfoil
![]()
图 12 不同网格量计算压力分布与试验压力分布对比
Figure 12. Comparison of pressure distributions among computational and experimental results
![]()
图 15 阻力系数代理模型预测误差收敛过程
Figure 15. The convergence history of surrogate prediction errors of drag coefficients
![]()
图 17 优化翼型压力分布对比
Figure 17. Comparison of pressure distributions of optimized and initial airfoils
![]()
图 18 单点优化设计翼型在非设计点气动特性对比
Figure 18. Comparison of aerodynamic characteristics of optimized airfoil on off-design points
![]()
图 20 优化翼型阻力发散特性对比
Figure 20. Comparison of drag-divergence characteristics of optimized and initial airfoils
![]()
图 21 优化翼型与RAE2822翼型在不同状态下压力分布对比
Figure 21. Comparison of pressure distributions of optimized and initial airfoils under different Mach numbers
![]()
图 23 不同计算网格量下计算压力分布与试验压力分布对比
Figure 23. Comparison of computational and experimental pressure distributions under different number of grids
![]()
图 25 两种方法阻力收敛历程对比
Figure 25. The convergence history of different optimization algorithms
![]()
图 26 初始机翼和优化后机翼(MF-PCK_VF-EI)上表面压力云图对比
Figure 26. Comparison of pressure distribution contours of upper surface of M6 wing using the MF-PCK_VF-EI method
![]()
图 27 初始机翼和优化后的机翼(MF-PCK_VF-EI)各剖面压力分布和翼型对比
Figure 27. Comparison of pressure distributions and shapes of initial and optimized wing sections
表 1 3种优化方法结果对比
Table 1 Comparison of optimization results using three methods
Optimization method The optimal
locationsThe optimal values Prediction errors/10−4 Kriging_EI [−3.1702579,12.2190345] 0.4174474 8.953 HK_VF-EI [3.1912351, 2.3199426] 0.4166618 7.642 MF-PCK_VF-EI [9.4185951, 2.4763138] 0.3981135 3.692 
下载: 导出CSV
表 2 RAE2822网格收敛结果对比
Table 2 Comparison of computational results of different levels of grids
Grid level Grid number $ {C}_{l} $ $ {C}_{d} $ $ {C}_{m} $ α/(°) $ {L}_{1} $ $4.0\times10^3$ 0.824 0.024297 −0.096164 3.022 $ {L}_{2} $ $ 2.1\times {10}^{4} $ 0.824 0.021261 −0.092284 2.973 $ {L}_{3} $ $ 8.1\times {10}^{4} $ 0.824 0.020112 −0.092868 2.917 $ {L}_{4} $ $1.27\times10^5$ 0.824 0.019960 −0.096062 2.869 $ {L}_{5} $ $1.98\times10^5$ 0.824 0.019729 −0.092778 2.806
下载: 导出CSV
表 3 RAE2822翼型优化设计结果对比
Table 3 Comparison of optimization results of RAE2822 airfoil
Optimization methods $ {C}_{l} $ ${C}_{d}/10^{-3}$ Areas ${C}_{m}/10^{-3}$ Total time of CFD evaluations Kriging_EI 0.824 11.375 0.0779 −91.98 $ 70{N}_{h} + 40{N}_{h} = 110 t $ HK_VF-EI 0.824 11.282 0.0779 −91.75 $ 50{N}_{h} + 15{N}_{h} + 500{N}_{l}\approx 90 t $ MF-PCK_VF-EI 0.824 11.182 0.0779 −91.71 $ 30{\mathit{N}}_{\mathit{h}} + 15{\mathit{N}}_{\mathit{h}} + 500{\mathit{N}}_{\mathit{l}}\approx 70\mathit{t} $
下载: 导出CSV
-
[1] Forrester AI, Keane AJ. Recent advances in surrogate-based optimization. Progress in Aerospace Sciences, 2009, 45(1-3): 50-79 doi: 10.1016/j.paerosci.2008.11.001
[2] Queipo NV, Haftka RT, Shyy W, et al. Surrogate-based analysis and optimization. Progress in Aerospace Sciences, 2005, 41(1): 1-28 doi: 10.1016/j.paerosci.2005.02.001
[3] Forrester DAIJ, Sóbester DA, Keane AJ. Engineering design via surrogate modelling: A practical guide. University of Southampton, 2008
[4] 韩忠华 许晨舟, 乔建领等. 基于代理模型的高效全局气动优化设计方法研究进展. 航空学报, 2020, 41(5): 623344 (Han Zhonghua, Xu Chenzhou, Qiao jianling, et al. Recent progress of efficient global aerodynamic shape optimization using surrogate-based approach. Acta Aeronautica et Astronautica Sinica, 2020, 41(5): 623344 (in Chinese) [5] Forrester AI, Sóbester A, Keane AJ. Multi-fidelity optimization via surrogate modelling. Proceedings of the Royal Society a:Mathematical, Physical and Engineering Sciences, 2007, 463(2088): 3251-3269 doi: 10.1098/rspa.2007.1900
[6] Peherstorfer B, Willcox K, Gunzburger M. Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Review, 2018, 60(3): 550-591 doi: 10.1137/16M1082469
[7] Palar PS, Tsuchiya T, Parks GT. Multi-fidelity non-intrusive polynomial chaos based on regression. Computer Methods in Applied Mechanics & Engineering, 2016, 305: 579-606
[8] Gano SE, Renaud JE, Sanders B. Hybrid variable fidelity optimization by using a kriging-based scaling function. AIAA Journal, 2005, 43(11): 2422-2433 doi: 10.2514/1.12466
[9] Kennedy MC, O'Hagan A. Predicting the output from a complex computer code when fast approximations are available. Biometrika, 2000, 87(1): 1-13 doi: 10.1093/biomet/87.1.1
[10] Han ZH, Görtz S. Hierarchical kriging model for variable-fidelity surrogate modeling. AIAA Journal, 2012, 50(9): 1885-1896 doi: 10.2514/1.J051354
[11] Keane AJ. Cokriging for robust design optimization. AIAA Journal, 2012, 50(11): 2351-2364 doi: 10.2514/1.J051391
[12] Toal DJ. Some considerations regarding the use of multi-fidelity Kriging in the construction of surrogate models. Structural and Multidisciplinary Optimization, 2015, 51(6): 1223-1245 doi: 10.1007/s00158-014-1209-5
[13] Zhao H, Gao ZH, Xu F, et al. Adaptive multi-fidelity sparse polynomial chaos-Kriging metamodeling for global approximation of aerodynamic data. Structural and Multidisciplinary Optimization, 2021, 64(2): 829-858 doi: 10.1007/s00158-021-02895-2
[14] Ng LWT, Eldred M. Multifidelity uncertainty quantification using non-intrusive polynomial chaos and stochastic collocation//Proceedings of the 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference 20th AIAA/ASME/AHS Adaptive Structures Conference 14th AIAA, 2012: 1852
[15] 赵欢, 高正红, 夏露. 基于新型多可信度代理模型的多目标优化方法研究. 航空学报, 2023, 44(2): 126962 (Zhao Huan, Gao Zhenghong, Xia Lu. Research on novel multi-fidelity surrogate model assisted many-objective optimization method. Acta Aeronautica et Astronautica Sinica, 2023, 44(2): 126962 (in Chinese) [16] Haftka RT, Villanueva D, Chaudhuri A. Parallel surrogate-assisted global optimization with expensive functions–A survey. Structural and Multidisciplinary Optimization, 2016, 54(1): 3-13 doi: 10.1007/s00158-016-1432-3
[17] Zhang Y, Han ZH, Zhang KS. Variable-fidelity expected improvement method for efficient global optimization of expensive functions. Structural and Multidisciplinary Optimization, 2018, 58(4): 1431-1451 doi: 10.1007/s00158-018-1971-x
[18] Huang D, Allen TT, Notz WI, et al. Sequential kriging optimization using multiple-fidelity evaluations. Structural and Multidisciplinary Optimization, 2006, 32(5): 369-382 doi: 10.1007/s00158-005-0587-0
[19] Jiang P, Cheng J, Zhou Q, et al. Variable-fidelity lower confidence bounding approach for engineering optimization problems with expensive simulations. AIAA Journal, 2019, 57(12): 5416-5430 doi: 10.2514/1.J058283
[20] 赵欢, 高正红, 夏露. 基于新型高维代理模型的气动外形设计方法研究. 航空学报, 2023, 44(1): 126924 (Zhao Huan, Gao Zhenghong, Xia Lu. Research on novel high-dimensional surrogate model-based aerodynamic shape de-sign optimization method. Acta Aeronautica et Astronautica Sinica, 2023, 44(1): 126924 (in Chinese) [21] Salehi S, Raisee M, Cervantes MJ, et al. An efficient multifidelity ℓ1-minimization method for sparse polynomial chaos. Computer Methods in Applied Mechanics & Engineering, 2018, 334: 183-207
[22] Schobi R, Sudret B, Wiart J. Polynomial-chaos-based Kriging. International Journal for Uncertainty Quantification, 2015, 5(2): 171-193 doi: 10.1615/Int.J.UncertaintyQuantification.2015012467
[23] Zhao H, Gao Z, Xu F, et al. An efficient adaptive forward–backward selection method for sparse polynomial chaos expansion. Computer Methods in Applied Mechanics & Engineering, 2019, 355: 456-491
[24] Jones DR, Schonlau M, Welch WJ. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 1998, 13(4): 455-492 doi: 10.1023/A:1008306431147
[25] 赵欢, 高正红, 夏露. 高速自然层流翼型高效气动稳健优化设计方法研究. 航空学报, 2021, 42(7): 124894 (Zhao Huan, Gao Zhenghong, Xia Lu. Research on efficient robust aerodynamic design optimization method of high-speed and high-lift NLF airfoil. Acta Aeronautica et Astronautica Sinica, 2021, 42(7): 124894 (in Chinese) doi: 10.7527/S1000-6893.2020.24894 [26] Zhao H, Gao Z, Xu F, et al. Review of robust aerodynamic design optimization for air vehicles. Archives of Computational Methods in Engineering, 2019, 26(3): 685-732 doi: 10.1007/s11831-018-9259-2
[27] 赵欢. 基于代理模型的高效气动优化与气动稳健设计方法研究. [博士论文]. 西安: 西北工业大学, 2020 Zhao Huan. Research on efficient surrogate-based aerodynamic optimization and robust aerodynamic design methods. [PhD Thesis]. Xi'an: Northwestern Polytechnical University, 2020 (in Chinese)
[28] Zhao H, Gao Z, Gao Y, et al. Effective robust design of high lift NLF airfoil under multi-parameter uncertainty. Aerospace Science and Technology, 2017, 68: 530-542 doi: 10.1016/j.ast.2017.06.009
[29] Zhao H, Gao Z. Uncertainty-based design optimization of NLF airfoil for high altitude long endurance unmanned air vehicles. Engineering Computations, 2019, 36(3): 971-996 doi: 10.1108/EC-05-2018-0215
[30] Zhao H, Gao Z, Xia L. Efficient aerodynamic analysis and optimization under uncertainty using multi-fidelity polynomial chaos-Kriging surrogate model. Computers & Fluids, 2022, 246: 105643
[31] 李焦赞 高正红. 基于几何不确定性的翼型多目标稳健优化设计. 力学学报, 2011, 43(3): 611-615 (Li Jiaozan, Gao Zhenghong. Multi objective optimization methodology for airfoil robust design under geometry uncertainty. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(3): 611-615 (in Chinese) doi: 10.6052/0459-1879-2011-3-lxxb2010-066 -
期刊类型引用(8)
1. 陈静，柯世堂，李文杰，朱庭瑞，员亦雯，任贺贺. 浅海风-浪-流-海床耦合场非定常时空演化规律及评价指标. 上海交通大学学报. 2023(06): 666-679 . 
百度学术
2. 李晓东，吴乘胜，王星. 船舶航行兴波CFD模拟的代数型VOF方法研究. 水动力学研究与进展A辑. 2022(02): 252-261 . 百度学术
3. 赵西增，徐天宇，谢玉林，吕超凡，姚炎明，解静，常江. 基于卷积神经网络的涵洞式直立堤波浪透射预测. 力学学报. 2021(02): 330-338 . 本站查看
4. 郭权势，邓争志，王晓亮，程鹏达. 垂荡双气室振荡水柱波能装置水动力特性研究. 力学学报. 2021(09): 2515-2527 . 本站查看
5. 吕超凡，赵西增，殷铭简. 波浪作用下涵洞式直立堤水体交换特性研究. 海洋工程. 2021(05): 39-49 . 百度学术
6. 柯世堂，王硕，张伟，王岩，王同光. 风、浪、流多场耦合作用波浪传播演化机理对比分析. 哈尔滨工程大学学报. 2021(09): 1312-1320 . 百度学术
7. 赵西增，聂隆锋，殷铭简. 一种高保真超长黏性数值波浪水槽. 水动力学研究与进展(A辑). 2020(01): 23-30 . 百度学术
8. 赵海洋，明平剑，张文平，齐文亮. 基于THINC/QQ格式的两相界面流动数值模拟. 哈尔滨工程大学学报. 2020(12): 1804-1810 . 百度学术
其他类型引用(3)
计量
- 文章访问数: 822
- HTML全文浏览量: 330
- PDF下载量: 154
- 被引次数: 11
