2. 洛桑联邦理工(EPFL), 瑞士 洛桑 CH-1015
在航空工程应用领域, 计算机数值模拟正在发挥着越来越重要的作用。但繁重的计算代价目前仍然制约着数值模拟的效率。为此, 在空气动力学、气动弹性等领域, 以杜克大学的Dowell、NASA的Silva等为代表的学者们广泛采用降阶模型(reduced-order model, ROM)[1-3]取代数值求解器来预测计算结果以减小计算量。
ROM是将高自由度全阶模型近似投影到子空间获得的低阶数学模型, 能以较少的自由度来描述全阶模型的主要动力学特性, 从而快速精确地给出全阶模型的预测解。目前应用最为广泛的ROM是通过本征正交分解方法(POD)提取流场的特征模态, 再通过适当的投影方法如Galerkin投影、一般投影等预测投影系数获得全阶流场的近似结果[4-5]。但是采用投影方法形式的降阶模型依赖于控制方程, 这给预测带来了不便, 如对于偏微分方程问题, 还需求解快照的梯度信息, 对于旨在减少计算量的ROM来说, 这种降阶模型组合显然不是最经济高效的。基于此, 研究学者们开始寻求建立系数求解不依赖于控制方程的非嵌入式降阶模型。邱亚松、白俊强等[6]利用POD和Kriging组成的非嵌入ROM成功预测了定常情况下的翼型流场。Fabien Casenave等[7]建立了基于离散经验插值方法(discrete empirical interpolation method, DEIM)的非嵌入式ROM, 并将其应用于气动噪声预测。而张伟伟等[8]则采用了POD与径向基函数法(Radial Basis Function, RBF)方法的结合预测翼型的气动响应和极限环振荡。
在ROM的实际应用中, 采用降阶模型可以预测出带参数问题在参数变化时的计算结果, 这对工程实践具有重要的意义。如邱亚松等[9]用一系列参数对翼型几何进行参数化, 然后通过降阶模型成功预测了不同翼型的定常流场。但目前降阶模型在带参数问题的应用大多是在定常条件下[10], 对于非定常的带参数问题则研究相对较少。考虑到同一个问题非定常计算量远高于定常, 此时降阶模型的收益应更可观, 另外对于某些精细化的设计如气动减噪等, 定常处理显得不够精确, 所以对带参数非定常问题建立降阶模型具有研究的价值与意义。Saifon Chaturantabut等[11]采用离散经验插值方法求解了一维Fitzhugh-Nagumo模型, 而N C Nguyen等分别对带参数非定常Boussinesq方程[12]和对流-扩散方程[13]建立了降阶基模型。但这些研究多对于简单方程, 对于空气动力学所研究的Navier-stokes方程则少有应用。
本文针对带参数的非定常偏微分系统, 通过采用基于本征正交分解方法(POD)的2层本征正交分解[14]方法提取带参数系统的特征模态, 再通过径向基函数(RBF)法预测各模态的投影系数, 建立了一种预测过程不依赖于系统所用控制方程的非嵌入式降阶模型。为检验降阶模型的精度与效率, 采用Burgers方程和不可压缩Naver-Stokes方程控制的驱动腔流动作为测试算例。
1 降阶基模型建立 1.1 模型建立对于本文研究的带参数的偏微分方程问题, 设目标函数uθ(x, t)x∈Ω, t∈[t0, t1], 边界条件和初始条件为:
(1) |
注意这里的参数θ是个矢量, 即该系统并不局限于单一参数。由于该问题涉及到空间域、时间域、参数域, 此时采用单一基函数φi(x)的传统降阶模型形式
(2) |
并不能直接处理, 本文引入两次正交分解(bi-orthogonal decomposition, BOD)[15], 其形式:
(3) |
式中,αθ(k, m)代表预测公式中待求系数, ϕk(x)、ξm(t)分别表征以空间坐标和时间作为自变量的降阶基, 简称为空间基函数、时间基函数。为便于构造, ϕk(x)、ξm(t)需要满足
(4) |
vθ(x, t)是为使预测结果
在求解空间、时间基函数ϕk(x)、ξm(t)即基模态时, 首先通过相关途径获取相应的快照(由Sirovich[16]首次提出):
(5) |
(6) |
快照的选取通常采用拉丁超立方[17]等随机采样方法。对这两组快照分别应用POD独立求解ϕk(x)、ξm(t)。这里简单介绍POD的基本方法。
对于一个非线性系统y(x, t), 系统提供N个快照yi(x)(x∈Ω, 1≤i≤N)。预测公式:
(7) |
Φi(x)(i=1, 2, …, N)表征待构造的单位正交基函数, 限定
为构造最小二乘意义上的最优基向量组, 快照的预测值
(8) |
‖·‖代表向量或者矩阵的二范数。通过转化, 可得到基函数Φi(x)所需满足的条件:
(9) |
即基函数应使快照在基函数张成的子空间内投影最大。定义快照的协方差矩阵M:
(10) |
要使快照投影达到最大, MΦi与Φi共线且同向, 因此只要求解矩阵M的正特征值, 便可以得到相应的基函数。而M是对称非负定矩阵, 故矩阵M具有N个非负特征值λi(i=1, 2, …, N):
(11) |
对求得的特征值λi(i=1, 2, …, N)进行排序, 再利用Sirovich[16]提出的广义能法选取前K个特征值λi及相应的特征向量vi。ε∈(0, 1)是给定的能量阈值, 表示选取的子空间蕴含的能量占总能量的比值, 通常给0.99。K的确定公式如下:
(12) |
由特征向量vi(x)可得到一组特征函数
(13) |
得到
采用POD方法, 利用快照Ψs、Ψt获得每个参数点θi的基函数ϕk, i(x)、ξm, i(t), 将所有参数点的基函数叠加, 形成模态集合, 最后应用SVD(奇异值分解)对基函数集进行二次压缩, 获得了最终的ϕk(x)、ξm(t)。此时, ϕk(x)、ξm(t)与θ无关, 可以用来预测参数域内任一点的目标函数。上述确定ϕk(x)、ξm(t)的整个过程被称为双层POD方法。
1.3 辅助函数求解为使结果满足边界条件和初始条件并保证降阶基的数学性质, BOD预测公式(3)中引入的辅助函数
(14) |
将辅助函数vθ(x, t)的求解分成v1θ(x, t)v2θ(x, t) 2个部分:
(15) |
(16) |
对这两个部分需要分别建立POD+Galerkin或者POD+RBF降阶模型进行预测, 由于这部分内容考虑的方程简单, 并且采用方法与之前有所重叠, 因此对预测过程不再进行详细介绍。但当边界条件不随时间变化时, 可以证明令vθ(x, t)=u0(x)时, 预测结果满足边界条件和初始条件:
当vθ(x, t)=u0(x)时
因此vθ(x, t)=u0(x)时,
要预测参数域内任一点对应的投影系数αθ(k, m)(1≤k≤K, 1≤m≤M), 首先需要利用已知的快照数据获得快照的投影系数αθi(k, m)(1≤k≤K, 1≤m≤M, 1≤i≤N)。为获得一组最接近样本点的投影系数集, 引入求最小值的优化问题:
(17) |
式中,‖·‖代表二范数, μ‖αθ(k, m)‖2项是为防止样本点的过拟合而添加的惩罚项。对于(17)式, 可通过带惩罚的最小二乘法[18]获得αθi(k, m):
(18) |
考虑到αθ是个以θ为自变量的函数, 每个元素是K×M的矩阵。并且通过最小二乘法的工作, 已有N组样本数据αθi(k, m), 对其采用RBF插值:
(19) |
φ(r)是RBF核函数, 本文采用高斯核函数:
(20) |
σ是高斯核函数的拟合参数, 需要人为给定。它对拟合结果影响非常大, σ绝对值过大, 影响拟合精度, 而σ绝对值过小, 又会容易造成拟合过程的失效。
令θ=θi(i=1, 2, …, N)及相应的αθi(k, m), 建立方程组求解γi(k, m)。
在求得γi(k, m)后, 对于任一个参数点θ, 直接通过(19)式, 便可算得其对应的投影系数αθ(k, m), 再由BOD公式(3)预测目标函数
至此, 降阶模型系统建立完毕。
2.3 RBF稳定性提高策略在实际应用过程中, 2.2节中展示的标准RBF在解决高维度问题时容易出现矩阵病态的情况, 并且对核函数中的拟合参数如高斯函数中的σ数值太过依赖。为此引入RBF-QR[19]来取代标准RBF。
RBF-QR并非是提出新的核函数φ(r), 而是通过了QR变换等数学手段对原有的核函数重构, 以尽可能地减小核函数中拟合参数σ对拟合效果的影响。
以本文采用的高斯核函数为例。通过泰勒展开, 坐标系转换和切比雪夫样条的引入, 将高斯核函数转变成:
(21) |
(22) |
式中, T(r)是切比雪夫样条组成的矩阵, D是对角矩阵, 矩阵元素为o(σ2j) (j=-1, -2, …, -∞)。矩阵C仅与核函数中自变量(r)相关。这部分过程主要是为了将矩阵中的r、σ相互分开。然后对矩阵C应用QR变换, 并构造新的核函数:
(23) |
DN、RN分别是矩阵D和R的前N×N个元素组成的矩阵, 拟合参数σ2j只存在矩阵D中, 在φ(r)的基础上乘以DN-1RN-1QT后, DN-1RN-1·R·D项的元素接近为o(1), 这样σ对新的核函数Ψ(r)的影响得到大幅度的减弱, 具体推导过程可见文献[19], 本文不作详细证明。
2.4 RBF-QR算例测试现采用2个简单的函数测试RBF-QR的性能。RBF核函数仍为高斯核函数。
a) 测试函数1
对本算例同时采用标准RBF(记为RBF-DIR)的拟合结果作为对比。取30个样本点βi=
2个自变量均在[0,1]范围内。测试时取200个样本点, 1 600个测试点。测试点的Halton分布如图 3a)所示, RBF-QR拟合的误差分布如图 3b)所示, 图 3说明在整个范围里, RBF-QR拟合结果都足够精确。
3 算例验证与结果分析本节对此双层POD+最小二乘法+RBF-QR的降阶模型进行算例验证。第一个算例是一维Burgers方程, 问题涉及到激波的预测; 第二个算例是不可压缩NS方程作为控制方程的二维驱动腔流动, 又考虑到了分离的捕捉。采用这两个具有典型非线性流动现象的算例对模型进行校核验证。
3.1 一维Burgers方程预测一维Burgers方程, 是模拟冲击波的传播和反射的非线性偏微分方程, 可用来研究激波的发展。方程形式如下:
(24) |
(25) |
υ=0.01, θ∈[-2,2]。这是个单参数并且边界条件不随时间变化的数学问题。采用有限元方法对偏微分方程迭代求解, 网格由100个均匀单元组成, 时间步长取为0.004 s, Nt=201。
在参数空间θ∈[-2,2]内均匀选取20个样本。令NT=51, Nx=51, 以此生成空间和时间快照。通过POD对每个样本点提取POD基模态后, 图 4展示了第10个样本点随模态数变化后提取的能量占总能量百分比的变化趋势, 很明显前几个模态占据了绝大多数能量。
令ε=0.001, 这样对所有样本点一共获得186个空间模态和207个时间模态, 再通过SVD方法对模态进行进一步压缩。再取ε=0.01, 最终得到K=8个空间模态和M=10个时间模态。
RBF-QR的参数取为σ=0.36。为评估模型的预测精度, 在θ∈[-2,2]内均匀选取N*=200个点作为测试样本。定义2个误差如下:
(26) |
降阶模型在这200个测试点的预测误差如图 5所示。
在整个参数范围内, 降阶模型的预测精度都比较高。图 6给出θ=1.94时4个时刻降阶模型和求解器得到的u分布的对比, 在对比图中2组曲线几乎完全吻合, 进一步证明了降阶模型在本算例中具有很高的预测精度。
3.2 驱动腔流动预测a) 问题描述
第二个验证算例为二维驱动腔流动, 驱动腔如图 7所示, 长度/m为1, 高度/m为h, 空腔上壁面施加均匀速度/(m·s-1)为θ, 下壁面及侧壁速度为0。θ和h作为本算例的2个可变参数, 给定它们的分布范围:0.5 m/s < θ < 2.5 m/s, 1 m < h < 2 m。
控制方程是不可压缩Navier-Stokes(NS)方程, 共有3个状态变量: x向速度U、y向速度V、压强p。二维不可压缩NS方程形式如下:
(27) |
在数值模拟时, 采用了一个采用压力修正的二阶有限差分求解器, 网格由51×51个均匀网格点组成, 时间步长0.01s, 计算由零初始条件(上边界除外)开始, 共预测T=5.0 s, 故Nt=501。为节省计算量, 生成空间快照时在[0, T]范围内均匀选取了NT=51个时刻, 而在生成时间快照的时候又用了Nx=11×11的均匀粗网格, 如表 1所示。
b) 降阶模型应用
考虑到这个问题仅仅涉及2个参数, 本文采用每个参数选取6个点的取样方式: h在1 m < h < 2 m内均匀选样; 而U则采用θ=
由于共有3个状态变量, 需要对它们分别提取降阶基。采用双层POD方法, 最终获得的模态数如表 2:
对3个状态变量分别预测时RBF-QR中的参数σ取为0.316, 0.316, 1。为评估降阶模型的性能, 选取3个测试点:(0.912, 1.3), (1.50, 1.7), (2.088, 1.7)(图 8中三角点)。定义降阶模型预测结果与真实值(数值模拟的结果视为真实值)的误差公式:
(28) |
在表 3中, 3个测试点的误差Error的绝对值都在0.01左右, 说明降阶模型对每个样本数据整体的预测精度高。但考虑到流场中仍有可能出现某些区域预测效果很差, 误差大的情况, 这也是我们不希望得到的结果, 因此观察测试点每个状态变量在整个区域内的误差分布, 如图 9所示。即使在预测最差的区域, 误差绝对值均小于0.1, 说明本算例中降阶模型对整个流场都实现了较精确的预测。
尽管误差分布的结果令人满意, 但在提取流线时, 本文发现用现有降阶基预测的流线分布与求解器得到的流线分布仍存在差距(如图 10a)、b))所示:样本二右下角大分离涡的上方存在一个小分离涡, 但降阶模型在此区域只能得到一个整体的分离涡; 对样本三右下角的大分离涡降阶模型捕捉不准, 没有模拟出涡核。为了准确捕捉流线分布, 将U、V的空间模态数均提高到30, 预测结果如图 10c)所示, 显然增加模态后, 降阶模型仿真出的流线分布与求解器的结果逼近程度高。
c) 计算效率评估
降阶模型的优势在于能够快速预测流场。本小节对它的整个计算效率进行了评估(如表 4所示, 注意这里降阶模型的耗时是增加模态数后的计算耗时)。
在得到样本点后更是降阶模型的训练阶段中提取每个样本点各自的基模态时共需要75.679 s, 对这些模态进行压缩时则需要3.844 s, 最后用最小二乘法求解αθi(k, m)需要27.65 s, 总耗时107.173 s。
模型建立后, 在预测某个点的流场结果时, 求解器需要130.42 s而降阶模型仅需要0.047 7 s, 计算耗时降为求解器的1/2734, 说明本降阶模型在驱动腔算例中不仅预测精度高, 而且非常高效。
4 结论与展望本文对带参数的非定常非线性系统建立了双层POD加最小二乘法加RBF-QR的非嵌入式降阶模型, 并得到以下结论:
1) RBF-QR相对于标准RBF方法不仅在相同σ下精度高, 而且更加稳定鲁棒, 对σ赋值要求也低;
2) 在一维Burgers方程和驱动腔流动中, 本文所采用的降阶模型能够处理非定常带参数问题;
3) 驱动腔内流算例中, 由于流线对结果非常敏感, 用采用能量法选取的降阶基并不能很好地复现流线分布, 这时仅增加速度的空间模态数就可以得到精确的流线。
4) 在驱动腔流动算例中, 预测一个样本, 降阶模型的计算效率相对于求解器的求解结果要提高2 734倍, 说明本文采用的降阶模型对研究非定常带参数问题具有非常可观的计算收益, 精确又高效。
但是, 为了避免自由度过多导致预测失效, 本文算例采用的是一维、二维参数空间的简单算例, 后续将考虑在保证降阶模型预测精度的前提下, 提高预测对象参数空间的维度和预测对象本身的复杂程度。
[1] | Dowell E H. Eigenmode Analysis in Unsteady Aerodynamics:Reduced Order Models[J]. AIAA Journal, 1996, 34(8): 1578-1583. DOI:10.2514/3.13274 |
[2] | Hesthaven J S, Rozza G, Stamm B. Certified Reduced Basis Methods for Parametrized Partial Differential Equations[M]. Springer, 2016 http://link.springer.com/978-3-319-22470-1 |
[3] | Quarteroni A, Rozza G. Reduced Order Methods for Modeling and Computational Reduction[M]. Springer, 2014 http://www.springerlink.com/openurl.asp?id=doi:10.1007/978-3-319-02090-7 |
[4] |
徐敏, 安效民, 曾宪昂, 等. 基于BPOD的气动弹性降阶及其在主动控制中的应用[J]. 中国科技论文, 2008, 10: 017.
Xu Ming, An Xiaomin, Zeng Xian'ang, et al. Model Reduction Using BPOD and Its Application to Aeroelastic Active Control[J]. China Sciencepaper, 2008, 10: 017. (in Chinese) |
[5] | Thomas J P, Dowell E H, Hall K C. Three-Dimensional Transonic Aeroelasticity Using Proper Orthogonal Decomposition Based Reduced Order Models[J]. Journal of Aircraft, 2003, 40(3): 544-551. DOI:10.2514/2.3128 |
[6] | Bai Jun Qiang, Qiu Ya Song, Qiao Lei. An Efficient Surrogate Model Construction Strategy for Large-Scale Output Problems[C]//2nd International Conference on Advances in Computational Modeling and Simulation, Kunming, China, 2013 |
[7] | Casenave F, Ern A, Lelivre T. A Nonintrusive Reduced Basis Method Applied to Aeroacoustic Simulations[J]. Advances in Computational Mathematics, 2015, 41(5): 961-986. DOI:10.1007/s10444-014-9365-0 |
[8] | Zhang W, Kou J, Wang Z. Nonlinear Aerodynamic Reduced-Order Model for Limit-Cycle Oscillation and Flutter[J]. AIAA Journal, 2016, 54(10): 3304-3311. DOI:10.2514/1.J054951 |
[9] | Qiu Yasong, Bai Junqiang. Stationary Flow Fields Prediction of Variable Physical Domain Based on Proper Orthogonal Decomposition and Kriging Surrogate Model[J]. Chinese Journal of Aeronautics, 2015, 28(1): 44-56. DOI:10.1016/j.cja.2014.12.017 |
[10] | Veroy K, Patera A T. Certified Real-Time Solution of the Parameterized Steady Incompressible Navier-Stokes Equations:Rigorous Reduced-Basis A Posteriori Error Bounds[J]. Int J Numer Methods Fluids, 2005, 47: 773-788. DOI:10.1002/(ISSN)1097-0363 |
[11] | Chaturantabut S, Sorensen D C. Discrete Empirical Interpolation for Nonlinear Model Reduction[C]//Proceedings of the 48th IEEE Conference on Decision and Control, 2009:4316-4321 |
[12] | Knezevic D J, Nguyen N C, Patera A T. Reduced Basis Approximation and a Posteriori Error Estimation for the Parametrized Unsteady Boussinesq Equations[J]. Mathematical Models and Methods in Applied Sciences, 2011, 21(7): 1415-1442. DOI:10.1142/S0218202511005441 |
[13] | Nguyen N C, Peraire J. An Efficient Reduced-Order Modeling Approach for Non-Linear Parametrized Partial Differential Equations[J]. International Journal for Numerical Methods in Engineering, 2008, 76(1): 27-55. DOI:10.1002/nme.v76:1 |
[14] | Christophe Audouze, Florian De Vuyst, Prasanth B Nair. Nonintrusive Reduced-Order Modeling of Parametrized Time-Dependent Partial Differential Equations[J]. Numerical Methods for Partial Differential Equations, 2013, 29(5): 1587-1628. DOI:10.1002/num.v29.5 |
[15] | De Wit T D, Pecquet A L, Vallet J C, et al. The Biorthogonal Decomposition as a Tool for Investigating Fluctuations in Plasmas[J]. Physics of Plasmas, 1994, 1(10): 3288-3300. DOI:10.1063/1.870481 |
[16] | Lawrence Sirovich. Turbulence and the Dynamics of Coherent Structures[J]. Quarterly of Applied Mathematics, 1987, 45(3): 561-571. DOI:10.1090/qam/1987-45-03 |
[17] | Sobol' I M. On the Systematic Search in a Hypercube[J]. SIAM Journal on Numerical Analysis, 1979, 16(5): 790-793. DOI:10.1137/0716058 |
[18] | Kaufman L. Maximum Likelihood, Least Squares, and Penalized Least Squares for PET[J]. IEEE Trans on Medical Imaging, 1993, 12(2): 200-214. DOI:10.1109/42.232249 |
[19] | Bengt Fornberg, Elisabeth Larsson, Natasha Flyer. Stable Computations with Gaussian Radial Basis Functions[J]. SIAM Journal on Scientific Computing, 2011, 33(2): 869-892. DOI:10.1137/09076756X |
2. Swiss Federal Institute of Technology in Lausanne(EPFL), CH-1015 Lausanne, Switzerland