張冬冬，郭勤濤

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

工程應用領(lǐng)域建模與仿真(M＆S)技術(shù)已被廣泛應用。NASA［1］頒布了建模與仿真的指導框架，旨在建立規(guī)范，促進其過程標準化。自建模與仿真技術(shù)產(chǎn)生以來，模型響應預測置信度一直為研究難點及熱點。國外對模型驗證與確認技術(shù)［2-4］進行研究，試圖應用于國防科研領(lǐng)域。ANS［5］制定出工程科學計算編程的模型驗證與確認指導框架;AIAA［6］發(fā)布計算流體力學仿真模型驗證與確認技術(shù)框架;ASME［7］發(fā)布計算固體力學仿真模型驗證與確認的技術(shù)框架。模型驗證與確認除含對軟件編程技術(shù)認證外，重點討論建模中各種不確定性［8-10］對模型響應不確定性影響?；谠囼炁c仿真數(shù)據(jù)確認方法［11-12］即如何判斷試驗與仿真預測間的一致性與相關(guān)性，也不斷探索與發(fā)展。盡管已有諸多研究，但模型驗證與確認因具體應用領(lǐng)域(或應用目的)不同，其運作思想存在諸多差別。為此，文獻［13-14］相繼提出模型驗證與確認的熱力學及結(jié)構(gòu)動力學挑戰(zhàn)問題，供全球?qū)＜覍W者探討模型驗證與確認的思路。此外，圣地亞國家實驗室正與美國能源部、國防部、工業(yè)界、學術(shù)界合作建立預測與狀態(tài)管理創(chuàng)優(yōu)中心，支持PHM技術(shù)的開發(fā)及技術(shù)試驗與確認。

吳立人［15］探討仿真模型有效性定量確認法，重點研究仿真模型有效性評估的定量化確認統(tǒng)計檢驗方法;郭勤濤等［16］初步探討了模型確認的主要內(nèi)容及總體框架，比較模型修正與模型確認的區(qū)別，并將模型確認技術(shù)應用于工程實例;鄧小剛等［17］總結(jié)計算流體動力學中模型確認方法及概念。目前，模型確認在國內(nèi)也逐漸受到重視，中國工程物理研究院與國家自然科學基金委員會［18］已于2008年批準“工程結(jié)構(gòu)數(shù)值模擬的模型驗證與確認方法研究”項目，重點研究工程結(jié)構(gòu)分層模型修正和確認技術(shù)框架，及模型驗證與確認的基本方法以及在工程中典型結(jié)構(gòu)中的應用。本文借鑒模型驗證與確認的基本思想，初步探討有限元模型確認流程，以Garteur Benchmark飛機結(jié)構(gòu)瞬態(tài)動力學仿真問題為例，采用三種不同響應面試驗設(shè)計方法，檢驗Kriging響應面的預測精度，借助Kriging響應面代理模型快速實現(xiàn)10萬個參數(shù)樣本的響應計算，采用核密度估計方法估計加速度響應最大值概率分布及置信區(qū)間上下限。

1 有限元模型確認流程

1.1 模型驗證與確認思想

模型驗證［7］定義為檢驗仿真計算模型的解是否與理論數(shù)學解準確吻合過程。模型驗證內(nèi)容包括:驗證仿真計算模型求解代碼編寫的正確性，算法實現(xiàn)的正確性，計算相對誤差及計算機代碼求解重復的可靠性。模型確認［7］定義為從目標用途角度出發(fā)，判斷仿真計算模型有多大可信度能準確描述真實世界。模型確認過程核心內(nèi)容是計算模型輸出響應結(jié)果與試驗響應結(jié)果的對比過程，據(jù)模型確認準則，實現(xiàn)計算模型可信度評價［19］。

圖1為模型驗證與確認流程圖。模型驗證與確認過程主要由兩個分支構(gòu)成，即建模分支與試驗分支。建模分支涵蓋三種模型:概念模型、數(shù)學模型、計算模型。例如，要預測某飛機機翼受外界氣動載荷作用影響產(chǎn)生的振動加速度，通常會將機翼簡化為梁，此簡化與假設(shè)過程即為建立概念模型過程。給定梁模型邊界條件，結(jié)合力學理論，建立梁振動求解方程，即為數(shù)學模型建立過程。計算模型通常以代碼形式出現(xiàn)，其基本特點是需將數(shù)學方程求解過程轉(zhuǎn)化為程序語言。模型驗證分為代碼驗證與精度驗證，代碼驗證通過對已知理論解問題的計算檢驗及減少計算模型在算法與編碼的誤差，精度驗證用于檢驗計算模型對物理試驗響應預測的準確度。試驗分支中設(shè)計的確認試驗為檢驗數(shù)學模型的計算精度與試驗不確定性的量化分析提供了數(shù)據(jù)基礎(chǔ)。

總之，模型驗證可保證計算模型代碼與算法的正確性、準確性、可靠性;模型確認可確認試驗響應與計算模型預測響應結(jié)果的量化對比，兩種數(shù)據(jù)結(jié)果的高度吻合性可表明概念模型假設(shè)簡化的可行性、數(shù)學模型與計算模型的正確性及準確性。

1.2 有限元模型確認流程

據(jù)模型驗證與模型確認基本思想，結(jié)合有限元方法的應用，本文提出有限元模型確認流程，如圖2所示。

圖1 模型驗證與確認流程圖［7］Fig.1 Flow chart of model verification and validation

圖2 有限元模型確認流程圖Fig.2 Flow chart of finite element model validation

有限元模型確認及技術(shù)問題的解決為:

(1)子結(jié)構(gòu)劃分及有限元模型建立

將復雜目標劃分為若干簡單子結(jié)構(gòu)，建立各子結(jié)構(gòu)有限元模型。據(jù)該模型參數(shù)是否存在不確定性，將其分為確定性子結(jié)構(gòu)有限元模型與不確定性子結(jié)構(gòu)有限元模型。

(2)校準試驗與確認試驗設(shè)計

校準試驗是通過試驗對子結(jié)構(gòu)模型參數(shù)，如材料彈性模量、泊松比、連接剛度等進行直接或間接測量及修正，統(tǒng)計分析不確定性參數(shù)。確認試驗是將記錄的試驗輸入數(shù)據(jù)及輸出響應數(shù)據(jù)作為不確定性子結(jié)構(gòu)有限元模型響應概率分布與試驗響應結(jié)果之間進行比較的依據(jù)。

(3)子結(jié)構(gòu)有限元模型確認

確定性子結(jié)構(gòu)有限元模型參數(shù)不確定性較小，可忽略不計。因此，該模型確認常采用確定性準則進行試驗與仿真響應結(jié)果直接比較。響應結(jié)果相對誤差較小時，接受該確定性子結(jié)構(gòu)有限元模型。而不確定性子結(jié)構(gòu)有限元模型參數(shù)不確定性較大，常采用不確定性準則，如置信區(qū)間準則，進行試驗與仿真響應結(jié)果比較;仿真樣本響應結(jié)果概率分布能較好涵蓋試驗響應值時，接受不確定性子結(jié)構(gòu)有限元模型，否則拒絕該模型。

(4)子結(jié)構(gòu)有限元模型合成

子結(jié)構(gòu)有限元模型合成技術(shù)方便了目標應用結(jié)構(gòu)有限元模型的建立，以該合成后模型作為目標應用結(jié)構(gòu)響應計算的有限元模型，可實現(xiàn)子結(jié)構(gòu)有限元模型參數(shù)向目標應用結(jié)構(gòu)有限元模型參數(shù)的直接轉(zhuǎn)化。

(5)代理模型建立

建立目標應用結(jié)構(gòu)有限元模型響應的代理模型，對其進行精度檢驗及偏差修正。工程中，單個目標應用結(jié)構(gòu)有限元模型計算時間長，因此，建立有限元模型的響應計算代理模型(或稱快速運行模型)。代理模型種類很多，如基于多項式的響應面代理模型、基于高斯過程的響應面代理模型、基于支持向量機的響應面代理模型等，其特點為可實現(xiàn)參數(shù)設(shè)計空間內(nèi)任意樣本點響應的快速、準確計算。

(6)參數(shù)蒙特卡洛抽樣計算

在目標應用結(jié)構(gòu)響應代理模型基礎(chǔ)上，結(jié)合蒙特卡洛抽樣計算方法，實現(xiàn)子結(jié)構(gòu)模型參數(shù)不確定性的正向傳遞。

(7)目標應用結(jié)構(gòu)響應預測

通常，目標應用結(jié)構(gòu)響應滿足一定概率分布，需采用統(tǒng)計學方法對其進行概率統(tǒng)計。典型的方法有參數(shù)型分布(如正態(tài)分布，指數(shù)分布)統(tǒng)計方法與非參數(shù)型分布統(tǒng)計方法(如核密度估計)。

2 Kriging響應面代理模型基本理論

代理模型(或稱元模型)在工程試驗設(shè)計中得到廣泛應用。為復雜結(jié)構(gòu)系統(tǒng)分析及優(yōu)化提供了方便。代理模型可描述為:給定n個參數(shù)樣本點并仿真試驗輸出，尋找能近似描述樣本點及仿真試驗輸出間關(guān)系的顯式數(shù)學模型。響應面模型作為代理模型具有操作簡單、能快速逼近等優(yōu)點。在計算機試驗設(shè)計分析(DACE)及優(yōu)化中，Kriging［20-22］插值方法起重要作用，以多項式逼近方式實現(xiàn)計算機仿真模型的近似描述，表達式［23］為:

其中:fj(x)為基函數(shù)，βj為基函數(shù)系數(shù)，z(x)為擬合偏差函數(shù)。Kriging方法認為不同樣本點處的擬合偏差量并非相互獨立，并假定該偏差函數(shù)為隨機過程Z(x)，且均值為0，方差為σ2，協(xié)方差非0。任意兩點t及u的協(xié)方差函數(shù)定義為:其中:ρ(t，u;θ)為相關(guān)函數(shù);θ為相關(guān)函數(shù)參數(shù)，用于衡量樣本點t與樣本點u之間相關(guān)性隨兩點間距離增加的衰減度，相關(guān)性參數(shù)越小，響應面越光滑。相關(guān)函數(shù)類型較多，通常以高斯相關(guān)函數(shù)作為常用相關(guān)函數(shù)。

其中:

式中:R為相關(guān)系數(shù)矩陣，Rij=ρ(xi，xj;θ)，(1≤i，j≤n)。

參數(shù)β，σ表達式為:

由此得:

相關(guān)系數(shù)θ采用數(shù)值求解方法，借助最大似然函數(shù)求極值確定其數(shù)值大小。記r(x)=［ρ(x，x1)，ρ(x，x2)，…，ρ(x，xN)］T，則 Kriging模型在任意設(shè)計參數(shù)樣點x處的響應預測表達式為:

3 實例應用

3.1 動力學問題及不確定性描述

圖3為Garteur benchmark飛機結(jié)構(gòu)示意圖。采用MD Nastran建立該飛機結(jié)構(gòu)有限元梁模型，機身與機翼，垂尾與平尾連接處均為固支，機翼兩端同時受幅值1 000 N、寬度1.1 ms的時域沖擊載荷F作用。設(shè)因加工制造工藝影響，飛機材料密度 ρ在區(qū)間［6.0，7.0］×10-9t/mm3滿足均勻隨機分布，彈性模量E在區(qū)間［2.0，3.0］×104MPa 滿足均勻隨機分布，采用該有限元模型預測測點在載荷作用方向加速度響應最大值概率分布。

圖3 Garteur benchmark飛機結(jié)構(gòu)示意圖Fig.3 The chart of Garteur benchmark structure

3.2 加速度響應最大值Kriging響應面

為檢驗響應面設(shè)計樣點數(shù)對Kriging響應面預測精度影響，選取拉丁超立方抽樣設(shè)計［24-25］、D最優(yōu)試驗設(shè)計、5水平全因子試驗設(shè)計及中心復合設(shè)計建立Kriging響應面模型見圖4，響應面檢驗點為5個隨機參數(shù)樣點，Kriging響應面精度檢驗結(jié)果見表1。比較四種試驗設(shè)計方法的 Kriging響應面精度檢驗指標R2、MSR、RMSE發(fā)現(xiàn)，拉丁超立方抽樣設(shè)計方法能以較少試驗點數(shù)實現(xiàn)設(shè)計空間內(nèi)響應量較準確預測，適用于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響應面以試驗設(shè)計參數(shù)樣本點為基礎(chǔ)，預測非試驗設(shè)計參數(shù)樣本組合下懸臂梁加速度幅值。結(jié)合蒙特卡洛及隨機抽樣方法，選8萬個彈性模量E與材料密度ρ的參數(shù)樣本組合，計算懸臂梁加速度幅值。加速度響應最大值非參數(shù)核密度估計［26］見圖5。由表2看出，在不同置信度條件下，加速度響應最大值的置信區(qū)間上下限并未隨置信度的下降發(fā)生顯著變化。

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

表2 不同置信度下置信區(qū)間Tab.2 Confidence interval under different confidence level

4 結(jié)論

通過討論模型驗證與確認，提出有限元模型確認流程;結(jié)合Kriging響應面理論，以Garteur benchmark飛機結(jié)構(gòu)有限元瞬態(tài)仿真為例，采用四種響應面設(shè)計方法建立加速度響應最大值Kriging響應面;借助蒙特卡洛和核密度估計方法的有效結(jié)合，實現(xiàn)飛機結(jié)構(gòu)加速度響應最大值概率分布預測;采用區(qū)間估計法，對飛機結(jié)構(gòu)加速度響應最大值置信區(qū)間上下限進行估計，結(jié)論如下:

(1)Kriging響應面代理模型能實現(xiàn)對有限元模型響應較準確預測，為有限元模型確認中模型參數(shù)不確定性傳遞，進一步實現(xiàn)目標應用結(jié)構(gòu)響應預測分析提供便利。

(2)置信區(qū)間準則能從統(tǒng)計學角度估計目標應用結(jié)構(gòu)響應上下限，當試驗響應數(shù)值大小位于置信區(qū)間上下限之間時，有限元計算模型可信度較好;否則，拒絕有限元計算模型。故可將置信區(qū)間準則作為有限元模型確認準則。

［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，F(xiàn)L，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］.導彈與航天運載技術(shù)，2000，4:24-26.WU Li-ren.Quantitative validation method of simulation model validity［J］.Missiles and Space Vehicles，2000，4:24-26.

［16］郭勤濤，張令彌.結(jié)構(gòu)動力學有限元模型確認方法研究［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］張保強，郭勤濤，陳國平.模型確認熱傳導挑戰(zhàn)問題求解的貝葉斯方法［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與響應面方法在氣動優(yōu)化設(shè)計中的應用［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模型在翼型氣動優(yōu)化設(shè)計中的應用研究［J］.西北工業(yè)大學學報，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模型的后優(yōu)化近似研究［J］.機械設(shè)計與研究，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］劉紀濤，劉 飛，張為華，等.基于拉丁超立方抽樣及響應面的結(jié)構(gòu)模糊分析［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統(tǒng)計分析與應用:40個案例分析［M］.北京:北京航空航天大學出版社，2010:355-369.