Kriging响应面代理模型在有限元模型确认中的应用

张冬冬，郭勤涛

(南京航空航天大学，南京 210016)

工程应用领域建模与仿真(M＆S)技术已被广泛应用。NASA［1］颁布了建模与仿真的指导框架，旨在建立规范，促进其过程标准化。自建模与仿真技术产生以来，模型响应预测置信度一直为研究难点及热点。国外对模型验证与确认技术［2-4］进行研究，试图应用于国防科研领域。ANS［5］制定出工程科学计算编程的模型验证与确认指导框架;AIAA［6］发布计算流体力学仿真模型验证与确认技术框架;ASME［7］发布计算固体力学仿真模型验证与确认的技术框架。模型验证与确认除含对软件编程技术认证外，重点讨论建模中各种不确定性［8-10］对模型响应不确定性影响。基于试验与仿真数据确认方法［11-12］即如何判断试验与仿真预测间的一致性与相关性，也不断探索与发展。尽管已有诸多研究，但模型验证与确认因具体应用领域(或应用目的)不同，其运作思想存在诸多差别。为此，文献［13-14］相继提出模型验证与确认的热力学及结构动力学挑战问题，供全球专家学者探讨模型验证与确认的思路。此外，圣地亚国家实验室正与美国能源部、国防部、工业界、学术界合作建立预测与状态管理创优中心，支持PHM技术的开发及技术试验与确认。

吴立人［15］探讨仿真模型有效性定量确认法，重点研究仿真模型有效性评估的定量化确认统计检验方法;郭勤涛等［16］初步探讨了模型确认的主要内容及总体框架，比较模型修正与模型确认的区别，并将模型确认技术应用于工程实例;邓小刚等［17］总结计算流体动力学中模型确认方法及概念。目前，模型确认在国内也逐渐受到重视，中国工程物理研究院与国家自然科学基金委员会［18］已于2008年批准“工程结构数值模拟的模型验证与确认方法研究”项目，重点研究工程结构分层模型修正和确认技术框架，及模型验证与确认的基本方法以及在工程中典型结构中的应用。本文借鉴模型验证与确认的基本思想，初步探讨有限元模型确认流程，以Garteur Benchmark飞机结构瞬态动力学仿真问题为例，采用三种不同响应面试验设计方法，检验Kriging响应面的预测精度，借助Kriging响应面代理模型快速实现10万个参数样本的响应计算，采用核密度估计方法估计加速度响应最大值概率分布及置信区间上下限。

1 有限元模型确认流程

1.1 模型验证与确认思想

模型验证［7］定义为检验仿真计算模型的解是否与理论数学解准确吻合过程。模型验证内容包括:验证仿真计算模型求解代码编写的正确性，算法实现的正确性，计算相对误差及计算机代码求解重复的可靠性。模型确认［7］定义为从目标用途角度出发，判断仿真计算模型有多大可信度能准确描述真实世界。模型确认过程核心内容是计算模型输出响应结果与试验响应结果的对比过程，据模型确认准则，实现计算模型可信度评价［19］。

图1为模型验证与确认流程图。模型验证与确认过程主要由两个分支构成，即建模分支与试验分支。建模分支涵盖三种模型:概念模型、数学模型、计算模型。例如，要预测某飞机机翼受外界气动载荷作用影响产生的振动加速度，通常会将机翼简化为梁，此简化与假设过程即为建立概念模型过程。给定梁模型边界条件，结合力学理论，建立梁振动求解方程，即为数学模型建立过程。计算模型通常以代码形式出现，其基本特点是需将数学方程求解过程转化为程序语言。模型验证分为代码验证与精度验证，代码验证通过对已知理论解问题的计算检验及减少计算模型在算法与编码的误差，精度验证用于检验计算模型对物理试验响应预测的准确度。试验分支中设计的确认试验为检验数学模型的计算精度与试验不确定性的量化分析提供了数据基础。

总之，模型验证可保证计算模型代码与算法的正确性、准确性、可靠性;模型确认可确认试验响应与计算模型预测响应结果的量化对比，两种数据结果的高度吻合性可表明概念模型假设简化的可行性、数学模型与计算模型的正确性及准确性。

1.2 有限元模型确认流程

据模型验证与模型确认基本思想，结合有限元方法的应用，本文提出有限元模型确认流程，如图2所示。

图1 模型验证与确认流程图［7］Fig.1 Flow chart of model verification and validation

图2 有限元模型确认流程图Fig.2 Flow chart of finite element model validation

有限元模型确认及技术问题的解决为:

(1)子结构划分及有限元模型建立

将复杂目标划分为若干简单子结构，建立各子结构有限元模型。据该模型参数是否存在不确定性，将其分为确定性子结构有限元模型与不确定性子结构有限元模型。

(2)校准试验与确认试验设计

校准试验是通过试验对子结构模型参数，如材料弹性模量、泊松比、连接刚度等进行直接或间接测量及修正，统计分析不确定性参数。确认试验是将记录的试验输入数据及输出响应数据作为不确定性子结构有限元模型响应概率分布与试验响应结果之间进行比较的依据。

(3)子结构有限元模型确认

确定性子结构有限元模型参数不确定性较小，可忽略不计。因此，该模型确认常采用确定性准则进行试验与仿真响应结果直接比较。响应结果相对误差较小时，接受该确定性子结构有限元模型。而不确定性子结构有限元模型参数不确定性较大，常采用不确定性准则，如置信区间准则，进行试验与仿真响应结果比较;仿真样本响应结果概率分布能较好涵盖试验响应值时，接受不确定性子结构有限元模型，否则拒绝该模型。

(4)子结构有限元模型合成

子结构有限元模型合成技术方便了目标应用结构有限元模型的建立，以该合成后模型作为目标应用结构响应计算的有限元模型，可实现子结构有限元模型参数向目标应用结构有限元模型参数的直接转化。

(5)代理模型建立

建立目标应用结构有限元模型响应的代理模型，对其进行精度检验及偏差修正。工程中，单个目标应用结构有限元模型计算时间长，因此，建立有限元模型的响应计算代理模型(或称快速运行模型)。代理模型种类很多，如基于多项式的响应面代理模型、基于高斯过程的响应面代理模型、基于支持向量机的响应面代理模型等，其特点为可实现参数设计空间内任意样本点响应的快速、准确计算。

(6)参数蒙特卡洛抽样计算

在目标应用结构响应代理模型基础上，结合蒙特卡洛抽样计算方法，实现子结构模型参数不确定性的正向传递。

(7)目标应用结构响应预测

通常，目标应用结构响应满足一定概率分布，需采用统计学方法对其进行概率统计。典型的方法有参数型分布(如正态分布，指数分布)统计方法与非参数型分布统计方法(如核密度估计)。

2 Kriging响应面代理模型基本理论

代理模型(或称元模型)在工程试验设计中得到广泛应用。为复杂结构系统分析及优化提供了方便。代理模型可描述为:给定n个参数样本点并仿真试验输出，寻找能近似描述样本点及仿真试验输出间关系的显式数学模型。响应面模型作为代理模型具有操作简单、能快速逼近等优点。在计算机试验设计分析(DACE)及优化中，Kriging［20-22］插值方法起重要作用，以多项式逼近方式实现计算机仿真模型的近似描述，表达式［23］为:

其中:fj(x)为基函数，βj为基函数系数，z(x)为拟合偏差函数。Kriging方法认为不同样本点处的拟合偏差量并非相互独立，并假定该偏差函数为随机过程Z(x)，且均值为0，方差为σ2，协方差非0。任意两点t及u的协方差函数定义为:其中:ρ(t，u;θ)为相关函数;θ为相关函数参数，用于衡量样本点t与样本点u之间相关性随两点间距离增加的衰减度，相关性参数越小，响应面越光滑。相关函数类型较多，通常以高斯相关函数作为常用相关函数。

其中:

式中:R为相关系数矩阵，Rij=ρ(xi，xj;θ)，(1≤i，j≤n)。

参数β，σ表达式为:

由此得:

相关系数θ采用数值求解方法，借助最大似然函数求极值确定其数值大小。记r(x)=［ρ(x，x1)，ρ(x，x2)，…，ρ(x，xN)］T，则 Kriging模型在任意设计参数样点x处的响应预测表达式为:

3 实例应用

3.1 动力学问题及不确定性描述

图3为Garteur benchmark飞机结构示意图。采用MD Nastran建立该飞机结构有限元梁模型，机身与机翼，垂尾与平尾连接处均为固支，机翼两端同时受幅值1 000 N、宽度1.1 ms的时域冲击载荷F作用。设因加工制造工艺影响，飞机材料密度 ρ在区间［6.0，7.0］×10-9t/mm3满足均匀随机分布，弹性模量E在区间［2.0，3.0］×104MPa 满足均匀随机分布，采用该有限元模型预测测点在载荷作用方向加速度响应最大值概率分布。

图3 Garteur benchmark飞机结构示意图Fig.3 The chart of Garteur benchmark structure

3.2 加速度响应最大值Kriging响应面

为检验响应面设计样点数对Kriging响应面预测精度影响，选取拉丁超立方抽样设计［24-25］、D最优试验设计、5水平全因子试验设计及中心复合设计建立Kriging响应面模型见图4，响应面检验点为5个随机参数样点，Kriging响应面精度检验结果见表1。比较四种试验设计方法的 Kriging响应面精度检验指标R2、MSR、RMSE发现，拉丁超立方抽样设计方法能以较少试验点数实现设计空间内响应量较准确预测，适用于Kriging响应面建立。

图4 加速度最大值Kriging响应面Fig.4 Kriging response surface for the maximum value of acceleration

表1 Kriging响应面精度检验Tab.1 Check for the accuracy of Kriging RSM

3.3 加速度最大值核密度估计

Kriging响应面以试验设计参数样本点为基础，预测非试验设计参数样本组合下悬臂梁加速度幅值。结合蒙特卡洛及随机抽样方法，选8万个弹性模量E与材料密度ρ的参数样本组合，计算悬臂梁加速度幅值。加速度响应最大值非参数核密度估计［26］见图5。由表2看出，在不同置信度条件下，加速度响应最大值的置信区间上下限并未随置信度的下降发生显著变化。

图5 加速度最大值核密度估计Fig.5 Kernel density estimation for the maximum value of acceleration

表2 不同置信度下置信区间Tab.2 Confidence interval under different confidence level

4 结论

通过讨论模型验证与确认，提出有限元模型确认流程;结合Kriging响应面理论，以Garteur benchmark飞机结构有限元瞬态仿真为例，采用四种响应面设计方法建立加速度响应最大值Kriging响应面;借助蒙特卡洛和核密度估计方法的有效结合，实现飞机结构加速度响应最大值概率分布预测;采用区间估计法，对飞机结构加速度响应最大值置信区间上下限进行估计，结论如下:

(1)Kriging响应面代理模型能实现对有限元模型响应较准确预测，为有限元模型确认中模型参数不确定性传递，进一步实现目标应用结构响应预测分析提供便利。

(2)置信区间准则能从统计学角度估计目标应用结构响应上下限，当试验响应数值大小位于置信区间上下限之间时，有限元计算模型可信度较好;否则，拒绝有限元计算模型。故可将置信区间准则作为有限元模型确认准则。

［1］NASA-STD-7009，National aeronautics and space administration［S］.Washington:NASA，2008.

［2］BalciO.Verification validation and accreditation for simulation models［C］.Proceedings of the 29th Conference on Winter Simulation，1997.

［3］ Aeschliman D P，Oberkampf W L.Experimental methodology for computational fluid dynamics code validation［R］.Sandia National Laboratories，Albuquerque，NM，1997.

［4］ Abgrall R，Desideri J A.The european hypersonic data base:a new CFD validation tool for the design of space vehicles［R］.AIAA 24th Fluid Dynamics Conference，Orlando，FL，1993.

［5］ ANSI/ANS-10.4-1987，American nuclear society:guidelines for the verification and validation of scientific and engineering computer programs for the nuclear industry［S］.

［6］G-077-1998，AIAA Guide for the verification and validation of computational fluid dynamics simulations［S］.

［7］ASME V＆V 10-2006，Guide for verification ＆ validation in computational solid mechanics［S］.

［8］McFarland J M.Uncertainty analysis for computer simulations through validation and calibration［D］.USA:Vanderbilt University，2008:136-156.

［9］Celik I，Zhang W M.Calculation of numerical uncertainty using richardson extrapolation:application to some simple turbulent flow calculations［J］.Journal of Fluids Engineering，1995，117(3):439-445.

［10］ Beck M B.Water quality modeling:a review of the analysis of uncertainty［J］.Water Resources Research，1987，23(8):1393-1442.

［11］ Bussoletti J E.CFD calibration and validation:the challenges of correlating computational model results with test data［C］.18th AIAA Aerospace Ground Testing Conference，Colorado Springs，CO，1994.

［12］ Sargent R G.Some approaches and paradigms for verifying and validating simulation models［J］.Proceedings of the Winter Simulation Conference，2001，1:106-114.

［13］ Dowding K J，Pilch M，Hills R G.Formulation of the thermal problem［J］.Computer Methods in Applied Mechanics and Engineering，2008，197(29-32):2385-2389.

［14］ Red-Horse R，Paez T L.Sandia national laboratories validation workshop:structural dynamics application［J］.Computer Methods in Applied Mechanics and Engineering，2008，197(29-32):2578-2584.

［15］吴立人.仿真模型有效性定量确认方法［J］.导弹与航天运载技术，2000，4:24-26.WU Li-ren.Quantitative validation method of simulation model validity［J］.Missiles and Space Vehicles，2000，4:24-26.

［16］郭勤涛，张令弥.结构动力学有限元模型确认方法研究［J］.应用力学学报，2005，22(4):575-578.

GUO Qin-tao，ZHANG Ling-mi.Finite elementmodel validation in structural dynamics［J］.Chinese Journal of Applied Mechanics，2005，22(4):575-578.

［17］邓小刚，宗文刚，张来平，等.计算流体力学中的验证与确认［J］.力学进展，2007，37(2):279-288.

DENG Xiao-gang，ZONG Wen-gang，ZHANG Lai-ping，et al.Verification and validation in computer fluid dynamics［J］.Advances in Mechanics，2007，37(2):279-288.

［18］张保强，郭勤涛，陈国平.模型确认热传导挑战问题求解的贝叶斯方法［J］.航空学报，2011，32(7):1202-1209.

ZHANG Bao-qiang，GUO Qin-tao，CHEN Guo-ping.Solution of modelvalidation thermalchallengeproblem using a bayesian method ［J］.Journal of Aeronautics，2011，32(7):1202-1209.

［19］ Rebbaan R，Mahadevan S.Computational methods for model reliability assessment［J］.Reliability Engineering and System Safety，2007，93(8):1197-1207.

［20］王晓峰，席 光，王尚锦.Kriging与响应面方法在气动优化设计中的应用［J］.工程热物理学报，2005，26(3):423-425.

WANG Xiao-feng，XI Guang，WANG Shang-jin.Application of kriging and response surface method in the aerodynamics optimization design ［J］. Journal of Engineering Thermophysics，2005，26(3):423-425.

［21］许瑞飞，宋文萍，韩忠华.改进Kriging模型在翼型气动优化设计中的应用研究［J］.西北工业大学学报，2010，28(4):503-509.

XU Rui-fei，SONG Wen-ping，HAN Zhong-hua.Application of improved kriging model based optimization method in airfoil aerodynamics design［J］. Journal of Northwestern Polytechnical University.2010，28(4):503-509.

［22］曹洪钧，段宝岩.基于Kriging模型的后优化近似研究［J］.机械设计与研究，2004，20(5):10-13.

CAO Hong-jun，DUAN Bao-yan.An approach on the post optimality approximation based on kriging model［J］.Machine Design and Research，2004，20(5):10-13.

［23］ Xiong Y.Using predictive models in engineering design:metamodeling，uncertainty quantification and model validation［M］.Evanston:Northwestern University，2008.

［24］吴振君，王水林，葛修润，等.LHS在边坡可靠度分析中的应用［J］.岩土力学，2010，31(4):1047-1054.

WU Zhen-jun， WANG Shui-lin， GE Xiu-run， etal.Application of latin hypercube sampling technique to slope reliability analysis［J］.Rock and Soil Mechanics，2010，31(4):1047-1054.

［25］刘纪涛，刘 飞，张为华，等.基于拉丁超立方抽样及响应面的结构模糊分析［J］.机械强度，2011，33(1):73-76.

LIU Ji-tao，LIU Fei，ZHANG Wei-hua，et al.Fuzz structure analysis based on latin hypercube sampling and response surface［J］.Journal of Mechanical Strength，2011，33(1):73-76.

［26］谢中华.MATLAB统计分析与应用:40个案例分析［M］.北京:北京航空航天大学出版社，2010:355-369.