李立
航空工业西安航空计算技术研究所,陕西 西安 710065
复杂航空工程取得成功的一个关键要素是在设计过程中引入数值计算(计算流体力学(CFD))的相关方法和工具。随着技术进步,CFD 已成为实现航空飞行器数字设计的关键性、基础性支撑工具。中国商飞在C919 飞机、C929飞机气动设计过程中均大量采用了CFD技术,通过将先进CFD、优化设计和试验验证等技术无缝结合,有力保障“设计具有较强竞争力的先进民用飞机”目标的实现[1]。
CFD的基本策略是综合利用计算机和数值算法来高效求解各种简化的或非简化的流体力学控制方程,获得关键性的空气动力学特性参数(如升力、阻力、载荷等)。从实际问题中抽象出来的流体力学控制方程是高度复杂的、高维非线性偏微分方程,目前对其数值解的数学理论研究并不充分。同时,由于实际物理问题本身的复杂性,CFD在具体应用时常常进行各种形式的数学简化。CFD技术的这些特殊性决定了针对具体问题,对数值模拟及结果的可信度进行综合分析和研判非常重要。目前,国际上公认的策略是对CFD 软件及其模拟结果进行验证(verification)、确认(validation)和不确定度量化(UQ)。
CFD 软件及模拟结果的验证、确认和不确定量化可以统一归属到CFD 可信度研究范畴。早在1998 年5 月,AIAA 杂志出版的可信的CFD 模拟(Credible CFD Simulations)中就明确提出,CFD可信度研究的几个关键问题包括验证、确认、认证和不确定度量化[2]。在美国国家航空航天局(NASA)和波音联合发布的、被视为当前CFD 技术发展风向标的CFD 2030愿景中,也明确把实现整机级问题的不确定度量化列为未来数十年CFD 重点研究的技术领域之一[3]。
本文在概述CFD 不确定量化方法当前发展现状的基础上,聚焦如何构建面向航空工程的CFD不确定度量化的高效算法及相应的支持平台和工具。具体涉及两个方面。一是方法层面的研究,包括标准术语、流程规范,以及瞄准旨在针对“维数灾难”,建立适于大维度不确定性输入参数不确定度量化的高效算法。二是支持平台和工具层面的研究,从开展第三方的验证、确认和不确定度量化思路出发,构建较为通用的不确定度量化基础平台框架。期望本文研究有助于助推国内相关工作开展。
在CFD 可信度分析范畴,不确定度量化是同验证、确认平行的基本概念,承担不同的职责和任务。1998年6月,基于不同行业专家联合研究成果,美国航空航天学会(AIAA)发布了CFD 验证和确认指南(AIAA G—077—1998)[4],被公认为是CFD可信度研究工作的权威指南。
该指南指出,建立CFD模拟可信度的两个基本原则是验证和确认。验证是指确定计算模型是否精确表示概念模型的过程,但对于计算模型或概念模型与真实世界的关系并无定论。确认是指确定计算模型和概念模型是否表示真实世界的过程。随后这一思路和概念持续得到发展。
2006 年和2009 年,美国机械工程师学会(ASME)先后发布计算固体力学验证和确认指南(ASME V&V10—2006)[5]、CFD 和热传输验证和确认标准(ASME V&V20—2009)[6]。这两个指南或标准的研究工作是近年在数值模拟软件及其模拟可信度研究相关术语、规范和标准研究中比较突出的工作。ASME V&V10—2006 指南和ASME V&V20—2009标准基本沿用了AIAA G-077—1998指南中其他相关术语和规范的定义,并继承了其软件可信度建立的思路和过程。
图1 给出ASME V&V20—2009 中提出的建立数值模拟软件可信度的基本过程。可以认为,目前国际上关于CFD 验证、确认和不确定度量化等基本术语及过程的认识趋于成熟。与AIAA G-077—1998 指南比较,明确提出除了验证和确认外,需将不确定度评估作为可信度分析过程的一项重要活动。
图1 建立数值模拟可信度的基本过程[6]Fig.1 Procedures to establish credibility for numerical simulation
2016 年,NASA 和波音在其联合发布的CFD 2030 愿景[4]中明确指出,需持续加强对CFD验证、确认和不确定度量化(VV&UQ)相关技术的研究,特别是在算法方面,到2018 年左右,在CFD 软件中应形成可靠的误差评估能力;到2025年左右,在CFD软件中应形成不确定度传播和对不确定度进行准确估计的能力;到2030 年,建立对所有可能来源的误差和不确定度(包括源于物理建模、数值误差、自然变异性和认知局限等方面的误差和不确定度)进行有效管理的能力,如图2所示。规划同时提出,在实际工作过程中,应关注不确定度量化技术与多学科分析和优化(MDAO)技术结合的问题。
图2 CFD 2030愿景中关于误差和不确定度量化的实施路径[4]Fig.2 Roadmap for error estimation and uncertainty quantification in CFD 2030 vision[4]
换言之,可以认为,实现CFD 全生命周期误差和不确定度的有效管理是对CFD 模拟进行验证和确认的终极目标。在开展CFD 模拟可信度分析过程中,引入“不确定度量化”这一工具,正是为了适应这种新趋势。
事实上,在复杂航空工程领域,CFD在实际使用过程中存在天然的不确定性。参考空中客车公司的思路和方法,从面向实际航空工程的角度出发,在CFD 模拟过程中,实际上主要存在三类典型的不确定性[7]:操作或运行不确定性、几何不确定性、物理建模不确定性。在复杂工程外形的计算评估或者优化设计问题中,如果建模时不考虑真实条件下上述输入条件的不确定性,往往无法准确评估复杂物理条件下真实外形的气动性能,尤其是考虑极限临界条件下的情况更是如此。在NASA 最新发布的报告中,不确定度量化被认为是支撑未来实现仿真认证的一项关键性技术[8],出发点即基于此。
针对不确定度量化相关课题,国外近10年开展了大量工作[9]。2007 年欧盟委员会在框架计划中设立了NODESIM-CFD 项目(基于CFD 非定值模拟的设计方法研究),随后继续在UMRIDA项目中进行推进;2007年12月,北大西洋公约下的RTO 实用飞行器技术小组专门组织内部专题技术会议,讨论军机设计中的不确定度量化问题;2009年,德国在其国家项目中设立了计算空气动力学中的误差和不确定度管理及优化(MUNA)项目,专门研究结合其主力计算软件Tau开展削减计算误差和不确定度的策略和方法等。
美国也设立有类似计划,通过各渠道项目支持不确定评估及管理相关工作开展。2013年Springer出版社编辑出版了名为“计算流体力学中的不确定度量化”的专门论著,报道了美国在CFD不确定度量化和管理方面的重要成果。主要包括三个方面内容:(1)CFD(包括气弹等非定常问题)中主要的不确定度评估问题;(2)高效、高精度的不确定度评估方法;(3)不确定度评估方法实现及验证。
2016年11月,在STI项目框架下,NASA组织有关团队对NASA近10年在验证、确认及不确定度管理方面开展的工作进行全面总结,编辑出版了《模拟的可信度:验证、确认和不确定度评估的研究进展》,基本体现了美国在相关研究的最高水平。主要内容包括:代码验证的策略和过程;计算验证的策略和过程;基于误差输运方程的解验证方法;自动化验证策略和方法;不确定度评估和敏感性分析的术语及规范;基于概率论的不确定度评估和敏感性分析方法;带随机参数的CFD计算及不确定度量化方法;试验设计方法及其应用;模型确认的方法;不确定度评估和验证、确认方法在试验及计算中的应用;高超声速吸气式推进系统计算确认的具体应用实践等。
在国内,笔者所在团队从2009年开始持续跟踪相关技术的研究,先后在国家科技部973计划课题(邓小刚院士任项目首席科学家)、装备预研重点基金、航空科学基金等渠道项目支持下开展工作,在基础理论和工具研发方面有一定经验和技术积累[10];北京应用物理研究中心在国家自然科学基金等项目支持下,针对核爆中的不确定度量化评估有关问题开展了研究,建立了针对核爆问题的一整套不确定评估办法[11];中国空气动力学研究与发展中心近年来在国家数值风洞工程支持下也启动了大量的工作,并以基础项目招标方式支持国内相关工作开展[12]。此外,北京航空航天大学、西北工业大学、国防科技大学等高校也结合相关任务,开展了CFD不确定度量化相关研究工作,限于篇幅,在此不再赘述。
应当说,目前,CFD不确定度量化问题逐渐引起国内相关团队的重视,但从研究的深度和广度,以及对不确定度量化问题本质的认识深入方面,笔者认为,与国外差距仍然相当明显,仍处于较初步的阶段。
不难分析,在开展CFD 可信度分析过程中,不确定度量化作为验证、确认下一阶段的工作,其主要工作目标是定量评估真实(数值)模拟条件下,不确定性因素对计算结果的影响。参考文献[13]对不确定度量化的基本内容及方法给出了很好的总结。
从工程实用化角度考虑,本文认为,CFD不确定度量化实际需要解决三个方面问题:(1)不确定性要素及误差来源的有效识别和定量刻画问题,即确定实际计算中有哪些不确定度要素来源会对实际模拟(或计算)产生大的影响;(2)不确定性如何传播、如何分析问题,即确定采用何种方法,实现各类不确定性在CFD 系统中的传播,并定量评估和确定不确定输入对计算结果有多大影响,以及在多大置信度区间范围内结果可信;(3)带不确定度的计算数据如何使用的问题,即有了含误差及不确定度区间的计算结果,如何基于不确定度量化结果进行设计决策。其中,最后一方面工作在实际研究中常常被忽略,但在工程上,这反而是最关心的。结合本文目的,下文重点论述前两个问题。
不确定性要素来源的有效识别和刻画是开展不确定度量化的第一个关键步骤,也是实现不确定度量化的前置条件。
在这一阶段,关心的主要问题包括:实际物理问题中,有哪些可能的不确定性影响因素?哪些因素对计算结果影响比较大,哪些因素在实际工作中可以忽略?也即确定哪些因素是关键因素、哪些因素是次要因素,从而实现在不确定性参数选取上,尽可能缩小参数的检索范围,减小发生“维数灾难”的可能性。然而,针对实际问题,上述不确定性要素来源的识别和刻画的工作本身并不总是特别容易完成。因而,从建立面向复杂航空工程的实用方法思路出发,围绕工程应用问题本身进行合理分析、建立恰当的不确定性要素分解图谱显得非常关键。
有大量文献针对CFD 不确定性影响因素进行如何分类开展研究[14]。在这些早期研究中,常把不确定度和误差混为一谈,将不确定性要素来源分为物理建模不确定性、应用不确定性、离散化和求解误差、计算机冗余错误、编程和用户错误等。但这种方法往往不利于分门别类建立恰当的不确定性传播与分析方法。主要原因是不确定性和误差在发生阶段、表征方式上均有相当大不同。可以认为,误差偏重于数值的部分,而不确定性偏重于物理的部分。
按照上述思路,在本文的具体实践中,拟从面向实际航空工程的角度出发,重点借鉴和参考空中客车公司曾采用的相关思路和方法,即面向复杂航空工程问题实际,重点关注下述三类不确定性要素,包括操作或运行不确定性、几何不确定性、物理建模不确定性[7]。
第一类是操作或运行不确定性,定义为因计算状态与实际运行状态有出入引起的不确定性。例如,与风洞试验对应,由于吹风条件或模型安装引起的流场均匀性、迎角等差异,往往会对模拟结果产生非常大的影响。按照传统分类方法,这是一类典型的应用不确定性要素。对于这类不确定度,一般认为,对于外流问题,马赫数、来流迎角和雷诺数被认为是最关键的影响因素,而对于内流问题,马赫数、来流迎角、环境压力和温度被认为是最关键的影响因素。
第二类是几何不确定性。定义为因几何模型制作工艺、模型表面粗糙度等导致的对真实几何外形描述有出入引起的不确定性。这也是一类典型的应用不确定度。对于这类不确定度,对于外流问题,模型的制作公差、模型的静态变形、模型的动态变形、温度影响(如结冰),以及持续的外在破坏因素(如雨滴影响)被认为是最关键的参数,而对于内流问题,模型的制作公差、模型的动态变形和持续的外在破坏因素被认为是最关键的参数。在具体计算中,这些参数的影响可以通过考察更具体的参数来加以评估。如对民机机翼模拟,值得考察的具体几何不确定度参数可能包括机翼前缘半径、机翼前缘和后缘后掠角、机翼尾缘厚度、机翼最大厚度及机翼最大厚度位置等。此外,气动弹性变形等也是值得关注的重要参数之一。
第三类是物理建模不确定性。定义为因采用的物理模型(包括所选用的几何模型和流动模型)与需要模拟的真实物理问题有出入引起的不确定性。按照这个定义,物理建模不确定度主要包含两个方面的问题。一类是几何建模的问题,这属于典型的应用不确定性问题。对于这类不确定性,几何模型是否简化被认为是最关键的参数。如在民机高升力机翼模拟中,为了便于网格生成和计算,往往会不具体模拟铰链和铆钉等部件,而采用干净的机翼模型。另外一类是流动建模的问题,这属于典型的模型不确定性问题。对于这类不确定性,明确各类流动模型的适用范围被认为是最关键的。对于湍流模拟,值得关注的具体关键参数是湍流模型中的模型常数和转捩参数。这些常参数往往通过试验结果标定,对不同具体流动,通常会带有很大的不确定性。
除上述三类不确定性要素外,数值参数引起的不确定性也是本文关注的。数值参数的不确定性可定义为因计算问题求解定义(包括初始条件、边界条件、空间离散格式、时间离散格式、计算网格生成以及相关的模拟参数)导致与实际有出入引起的不确定性。这部分问题既包含部分应用不确定性问题,如初始条件、边界条件的选择和设置,以及计算条件的物理属性参数设置等,又包含部分数值离散误差或求解误差引起的问题,如计算网格生成、空间离散格式和时间离散格式选择,以及相关求解参数设置等。
按照上述方法,实际上基本建立了一套适于航空工程的不确定性要素分解的图谱。
在不确定性的具体表征上,一般方法是按不确定性的具体表现形式分为随机不确定性和认知不确定性,进而分类表示[13-14]。随机不确定性是一类客观存在的不确定性,常采用概率密度函数(PDF)、均值、方差、偏度等概率统计理论来表达。认知不确定性是一类主观的不确定性,主要是由于对问题的认识不足造成的,因而通常很难像随机不确定性表征为概率密度函数,常常用非概率的方法来表征,如可选值、区间等。
不确定性传播与分析可以看作利用CFD 系统进行多参数学习、收集分析证据的过程。因此,不确定性传播与分析方法在不确定度量化中起到关键性作用。
依据不确定性表现形式不同,主要包括适用于随机不确定性的概率类方法(如蒙特卡罗(MC)方法、多项式混沌(PCE)方法、矩方法),以及适于认知不确定性的非概率类方法(如区间分析方法、模糊逻辑方法)[13]。此外,以伴随方程灵敏度导数方法为代表的敏感性导数方法也常用于认知不确定性传播与分析的典型方法[15]。
这里重点介绍两类方法。一类是以蒙特卡罗(MC)方法及其改进为代表的随机采样的方法;另一类是以非介入式多项式混沌(NI-PCE)为代表的结构性采样方法。这两类方法是笔者所在课题组建立大维度可变输入复杂问题模拟的高效不确定评估方法的基础。具体目标是,针对经典方法从低维拓展到高维带来的采样点数量呈几何级数增加引起的所谓“维数灾难”问题,开展方法的优化及改进。
2.2.1 经典蒙特卡罗及其改进方法
按照数学的观点,随机不确定性的传播和分析本质上可以看作一个随机微分方程的求解问题,因而,利用经典蒙特卡罗方法(MC)是一种非常直接、自然的思路。并且,在理论上,利用蒙特卡罗方法总能无限逼近随机微分方程的真解。从这个意义上来说,经典蒙特卡罗(MC)方法可以看作一种能对随机不确定性量化评估结果正确性进行验证的可靠方法。
MC方法的基本思路是:针对拟求解的随机问题,建立随机过程,通过对该随机过程进行采样和数值试验,得到拟求解问题的近似解。MC方法的主要优点是无须显式推导或给出影响变量与所研究变量之间的关系,而是基于对实际仿真数据的统计分析,得到结果的统计特征,具有简洁、明了的特点。然而,在实际应用中,MC方法最主要的缺点是收敛非常慢,其估计误差的精度可表示为
式中:N为采样样本数,δ为利用随机样本输出计算的方差,λα是与置信度α一一对应的常数,可以通过查表法得到。
由式(1)不难分析,传统MC 方法的收敛速度与N成反比:当N= 100时,模拟误差精度为10-1量级,要达到10-2量级,样本数至少增加到10000,效率非常低。因此,针对复杂工程问题,如何改进MC方法、提高其效率一直是重要研究之一。
在具体实践中,利用MC方法基于随机抽样的特点,一种有效的策略是根据变量特点减小随机抽样样本点数量。拉丁超立方蒙特卡罗方法(LHS-MC)就是典型的方法之一。但是该类型方法尽管在一定程度上可以提高计算效率,但并没有改变抽样类方法的本质。为了提高评估精度,仍需要大量计算。
在本文课题组的实践中,提出一种基于飞行仿真中数据融合的思路对MC方法进行改进的算法[16]。其主要思路是在仿真阶段,通过小样本计算,引入Kriging 代理模型来重构产生大样本,实现对MC 方法的改进,由此得到基于Kriging 代理模型的改进MC 方法,简称Kriging-MC 方法。步骤如下:(1)利用试验设计方法或伪随机方法产生给定分布的大样本初始输入;(2)利用均匀随机采样、拉丁超立方采样等抽样方法从该初始样本随机抽取给定样本容量的小样本作为真实计算输入;(3)进行重复性数值试验,得到小样本的计算结果;(4)基于上述小样本计算结果,利用Kriging响应面方法进行插值,得到与初始大样本随机输入对应的计算结果;(5)基于上述计算结果进行统计分析,给出具体的不确定度量化结果。
在上述过程中,采用Kriging响应面方法进行大容量样本数据的恢复方法,相较多项式响应面方法具有明显优势。Kriging 方法作为一类统计学意义最优的方法,已被广泛应用到数据融合、优化设计等领域。具体理论公式及分析可参见参考文献[16]。
2.2.2 非介入式多项式混沌方法
多项式混沌(PC)方法是近年来发展非常迅速的方法。PC 方法本质是一种谱分解方法,基于结构性采样。PC 方法将整个数值模拟过程认定为一个随机过程,根据输入参数的概率密度分布函数,选择合适的正交基函数,然后将该随机过程的输出在基函数构成的谱空间中展开,进而将问题核心转换为谱空间中正交多项式的系数确定问题。
PC方法从原理上就与MC方法有本质不同,由于利用谱空间的正交分解,计算量与期望的多项式的展开阶数正相关,因而计算效率非常高。
PC方法基本原理概述如下。
假定随机的输入变量ξ=(ξ1,ξ2,…,ξn)满足某种特定概率密度分布f(ξ)。对于任意关心的输出量y=y(ξ)(根据具体情况,可以是升力、阻力等总体积分量,也可以是流场中某处的压力系数、摩擦力系数或者马赫数等流场局部变量,或者分离区长度、再附点位置等流场特征量),利用谱分解理论,可以将其进行正交展开
式中:αj为待确定的系数;Ψj是以f(ξ)为权系数、具有随机性质的正交多项式,如对多维正态分布输入,对应混沌正交多项式为Hermite多项式;NPC是进行谱分解截断的总项数,可以依据混沌多项式的阶数p及随机变量个数NPC来确定
这样,问题转换为如何确定待定系数αj。一旦系数确定,则可以利用其计算任意输入量,并推导出输出量的统计特征
上述过程依据其与求解器耦合实现的方式,进一步可以分为介入式多项式混沌(IPC)方法及非介入式多项式混沌(NIPC)方法。
在IPC中,控制方程中的不确定变量和参数都用其PC展开替代,经过投影之后,得到与原控制方程相似的随机控制方程,自变量为PC展开中的自由度。IPC属于一种白盒方法,需要对现有代码做大的改动,基本上可视为开发了一种新的求解器,不仅理论推导难度大,而且程序改造工作量较大,加上很多情况下,软件代码都是封闭的,这限制了IPC在工程问题中的应用。
因此,从面向航空工程实际应用角度出发,本文主要发展和采用NIPC 方法。NIPC 方法是一种典型的黑盒方法,将CFD 求解器作为“黑箱”,通过确定性样本,对展开系数进行估计。确定展开系数的过程,也常称为匹配过程。主要有两种思路。一种称为Galerkin匹配法,依据式(2)可得
式(2)可以采用高斯积分方法进行计算。另一种称为配置法,采用任意多个样本的确定性计算结果进行待定系数确定。通常情况下,计算样本数N应大于PC展开的总项数NPC。这样,利用N个观测值yi(i= 1,2,…,N)可以形成一个关于待确定系数的超定方程:
求解上述方程即得到待确定的系数αj。
从发展现状来说,尽管与MC 方法相比,PC 方法的计算效率通常表现为上百倍效率的提升,但针对实际工程问题,随着维数增加,PC 方法同样存在展开的项数随着输入变量的维数和展开阶数的增加成几何级数增长、引起“维数灾难”的问题。其主要原因是在PC 方法中,通常采用张量积方式实现多维求积。
为此,人们在PC方法基础上,近年来开展了诸多尝试,期望进一步改进和提升方法的效率,以应用到三维复杂工程问题。稀疏网格重构多项式混沌(Sparse-PC)是其中值得关注的方法[17]。在包含多输入变量的随机问题中,经常会遇到目标响应的展开式稀疏的情况,展开中占主导的通常是独立的输入变量和变量之间的低阶交叉项。这样,某些非常小的自由度的量值完全可以舍去,这使得稀疏多项式混沌方法得以广泛应用。研究表明,通过引入稀疏网格技术,可以极大程度减小高维问题计算所需节点数,因而可在一定程度上缓解“维数灾难”问题。
为了面向航空工程应用,形成实战化能力,自主开发和研制具有较强通用性及平台适应能力的不确定度量化工具非常必要。目前,国际上已发布多种开源、可支持不确定度量化分析通用工具或平台。但总体针对性不强。为此,本文提出建立一种针对性及第三方软件适应性较强的不确定度量化评估基础平台框架UQ-CFD V1.0.
UQ-CFD V1.0采用C++/Fortran语言混合开发。其中,C++用来实现功能接口的封装、调用,以及与求解器通信,实现具体UQ 功能;Fortran 用来实现求解器前后处理的标准化,实现CFD 计算数据(力系数数据、表面场数据、空间场数据)的提取,方便与C++混合编译,为前端UQ分析程序自动提供监控数据。当前版本已实现并支持的主要UQ功能包括:随机样本自动生成(正态分布、均匀分布等)、基于Kriging 代理模型的改进MC 不确定度评估功能、基于稀疏网格积分的PC不确定度评估功能等。
图3 给出了UQ-CFD 的系统组成。包括UQ-PRE、UQ-UPDATE、UQ-POST 三个核心模块,分别承担不同任务。UQ-PRE 负责前处理,UQ-UPDATE 负责计算评估和结果更新,UQ-POST负责后处理。UQ-CFD致力于建立一个通用性的交互式不确定度量化评估应用环境及平台框架,提供单参数和多参数随机计算,并实现对LHS-MC、Kriging-MC、PC、Sparse-PC等高效率不确定度传播方法的支持,为开展复杂航空工程问题的CFD不确度量化奠定了基础。
图3 UQ-CFD框架的实现技术思路Fig.3 Technical roadmap for UQ-CFD implementation
为提高该工具的使用性,UQ-CFD提供了一种基于标准化参数集成思路的求解器耦合方法,通过引入一种较通用的参数标准,便捷地为实现不同求解器耦合带来的参数传递问题。通过与多个内部求解器的耦合验证,证实该策略非常有效。UQ-CFD提供两种运行模式:命令行交互模式及批处理脚本模式。图4展示了平台运行界面和主要系统功能。提供的主要系统功能包括:提供交互方式UQ应用环境;支持基于MC、PC 的典型UQ 方法;支持单变量和多变量随机样本生成;支持点数据、表面流场数据、空间全流场数据自动收集;支持自动运行脚本(C语言实现)等。
图4 UQ-CFD框架运行界面及主要系统功能Fig.4 System function and usage method of UQ-CFD
为了对支持的基础算法进行验证,采用经典Rosen brock函数作为测试案例。该函数定义为
确定性问题为:求函数f(ξ1,ξ2)在(0,0)处的值。对应的不确定度量化问题是:当ξ1,ξ2均为符合正态分布N(0,1)的随机输入时,求f(0,0)在95%置信度下的不确定度。
对该二维问题,采用5种方法进行计算:(1)标准MC方法;(2)LHS-MC 方法;(3)Kriging-MC 方法;(4)Galerkin-PC 方法;(5)Sparse-PC 方法。计算中,样本数统一设定为1000,对MC 类方法,为直接样本数;对PC 类方法,为重构样本采样数。为了便于对结果正确性进行检验,同时采用LHS-MC 进行了10000 次样本的大容量计算,将其视作基准结果。详见表1。
表1 不同方法得到的Rosenbrock函数在(0,0)处的统计特征Table 1 Comparison between statistics for Rosenbrock test function with different UQ approaches at(0,0)
此外,Kriging-MC 依据重构样本数不同,选择三种具体方法,以评估重构样本数对不确定度评估结果的影响;Galerkin-PC 依据阶数不同,选择两种方法,以评估阶数影响;Sparse-PC 依据层数选择两种方法,以评估稀疏网格层数影响。
对本算例,以均值、方差、偏度为主要判断标准,不难分析:(1)标准MC及LHS-MC方法给出的计算结果相当,结果合理;对本算例,LHS-MC 方法给出的偏度偏小。(2)基于Kriging代理模型的改进蒙特卡罗方法,重构样本数的增加带来预测准确性的提升,但同时带来计算量的显著增加(与重构样本数正相关,重构样本数20,需要计算20次);对本算例,重构样本数20,与重构样本数100、500 的方法,给出的结果已大致相当,表明方法的高效性。(3)对Galerkin-PC,阶数提升带来预测准确性的提升,尤其是高阶的偏度估计,但计算量相应增大(对本算例,5 阶方法需要25 次计算、10阶方法需要100次计算);以均值、方差为基准,5阶方法与LHS-MC 方法计算结果相当。(4)对Sparse-PC,层数增加同样会带来准确性提升,对本算例,2层Sparse-PCE已经能够给出相当准确的结果,3 层Sparse-PC 则与10 阶Galerkin-PC 给出的计算结果相当,但计算量只有10 阶Galerkin-PCE的一半(直接计算45次;10阶Galerkin-PC计算100次),因而计算效率有显著提升。
上述测试结果,证实了本文所采用基础不确定度量化方法的可靠性。
来流马赫数3.0 的平面斜激波问题是参考文献[15]中用于几何不确定性评估的经典案例。初始平面倾角为5o,在超声速来流作用下将产生斜激波,激波位置对平面顷角非常敏感。
对该问题,采用PC 方法进行不确定度量化分析,并采用随机样本1000 的标准MC 方法作为基准结果进行对比。图5给出来流马赫数3.0条件下,平面斜激波几何影响不确定度评估的典型结果,清晰地体现出几何(倾角5o±1%)对该流场的影响。由计算结果不难分析,该结果与理论分析表现出的特征完全一致。
图5 平面斜激波问题几何不确定度量化典型结果Fig.5 Typical UQ results for geometry uncertainty qualification of oblique shock problem
通过图6 可以看出,随着混沌多项式阶数升高,PC 与MC 方法两者的相对差量越来越小。对该问题,当混沌多项式的阶数达到6 时,两者评估结果已达到基本一致,表明PC方法相对MC方法的高效性。
图6 平面斜激波问题中,多项式阶数对PC评估结果的影响Fig.6 Effect of polynomial order on PC for oblique shock problem
本节重点在于考察多参数情况下的综合不确定度评估。确定性算例为经典的高超声速圆柱算例,计算条件为:马赫数8.03,自由来流温度124.94K,壁温294.44K,雷诺数为1.835e+5。
对应的不确定度量化问题定义为:来流马赫数、自由来流温度,以及壁温三个环境参数对模拟结果的影响。各不确定性参数统一采用1%不确定性,均为正态分布。
本算例选择Galerkin-PC 作为主要的不确定度量化方法。混沌多项式阶数取4,直接重复计算的次数为64次,与蒙特卡罗方法动辄需要数千次计算相比,计算量显著下降。图7对比了考虑不确定性计算得到的压力平均场与直接计算基准场,表明采用Galerkin-PC 给出的平均场结果合理,可作为进一步分析的基础。
图7 不同计算得到的压力场对比Fig.7 Comparsion between pressure contour result with different approaches
图8进一步给出对压力场标准差、偏度、斜度等统计分析情况。其中,标准差可以认为是流场影响量的最佳度量,其大小可以度量预测结果的可信度水平。可以看到,以标准差预测作为评判标准,对本算例,操作不确定性主要集中体现为对弓形激波位置的模拟可信度影响较大。
图8 考虑来流不确定性的高超声速圆柱问题压力场标准差、偏度、斜度结果Fig.8 Results for stand deviation,skewness,and kurtosis for hypersonic cylinder considering operational uncertainty
图9 进一步给出表面热流的计算情况。可以看到,按照95%置信度分析,考虑不确定性影响的计算方法很好地给出了热流计算结果。最大热流与试验吻合较好。
图9 考虑来流不确定性的高超声速圆柱热流结果Fig.9 Heat flux result for hypersonic cylinder considering operational uncertainty
利用上述结果,可以进一步采用Sobol灵敏度指数分析方法,对各参数影响量进一步定量分析。对本算例,灵敏度分析结果表明(见表2),以壁面阻力系数作为考察对象,计算结果受壁温影响最大(贡献度为49.9%)、来流静温次之(贡献度为31.6%),马赫数贡献最小(贡献度为18.5%)。
表2 高超声速圆柱问题灵敏度分析结果Table 2 Sensitivity analysis results for hypersonic cylinder problem
瞄准建立实用的方法和工具,本文对面向航空工程的CFD不确定度量化方法进行论述。本文研究得到以下几个结论:
(1)CFD 不确定度量化是可信度分析范畴中与验证、确认并驾齐驱的重要技术领域,随着人们对CFD技术需求的提升,不确定度量化在支撑鲁棒性数字设计及设计认证方面大有作为。
(2)瞄准工程实用性,建立适于大变量可变输入的不确定性量化方法非常关键,基于谱展开理论的多项式混沌(PC)及其优化及改进方法是最值得重点关注的方法。
(3)从构建实战能力出发,自主开发和研制具有较强通用性及平台适应能力的不确定度量化工具也非常必要,本文展示的UQ-CFD 提供了一种通用性不确定度量化工具平台框架的示范。
(4)基于标准Rosenbrock函数算例及两个典型应用案例的计算,表明本文的方法和策略是有效的,可以有效应用与几何不确定性、操作不确定性、物理模型不确定性等随机不确定性影响评估相关的问题。
目前,本文工作仍处于较初步阶段,下一阶段将进一步对研制的UQ-CFD框架功能进行完善,重点聚焦建立针对同时考虑随机不确性及认知不确定性的混合不确定度量化评估的能力,以及适应复杂问题的能力。