張冬冬，郭勤濤

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

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

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

1 有限元模型確認(rèn)流程

1.1 模型驗證與確認(rèn)思想

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

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

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

1.2 有限元模型確認(rèn)流程

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

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

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

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

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

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

(2)校準(zhǔn)試驗與確認(rèn)試驗設(shè)計

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

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

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

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

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

(5)代理模型建立

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

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

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

(7)目標(biāo)應(yīng)用結(jié)構(gòu)響應(yīng)預(yù)測

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

2 Kriging響應(yīng)面代理模型基本理論

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

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

其中:

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

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

由此得:

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

3 實例應(yīng)用

3.1 動力學(xué)問題及不確定性描述

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

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

3.2 加速度響應(yīng)最大值Kriging響應(yīng)面

為檢驗響應(yīng)面設(shè)計樣點(diǎn)數(shù)對Kriging響應(yīng)面預(yù)測精度影響，選取拉丁超立方抽樣設(shè)計［24-25］、D最優(yōu)試驗設(shè)計、5水平全因子試驗設(shè)計及中心復(fù)合設(shè)計建立Kriging響應(yīng)面模型見圖4，響應(yīng)面檢驗點(diǎn)為5個隨機(jī)參數(shù)樣點(diǎn)，Kriging響應(yīng)面精度檢驗結(jié)果見表1。比較四種試驗設(shè)計方法的 Kriging響應(yīng)面精度檢驗指標(biāo)R2、MSR、RMSE發(fā)現(xiàn)，拉丁超立方抽樣設(shè)計方法能以較少試驗點(diǎn)數(shù)實現(xiàn)設(shè)計空間內(nèi)響應(yīng)量較準(zhǔn)確預(yù)測，適用于Kriging響應(yīng)面建立。

圖4 加速度最大值Kriging響應(yīng)面Fig.4 Kriging response surface for the maximum value of acceleration

表1 Kriging響應(yīng)面精度檢驗Tab.1 Check for the accuracy of Kriging RSM

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

Kriging響應(yīng)面以試驗設(shè)計參數(shù)樣本點(diǎn)為基礎(chǔ)，預(yù)測非試驗設(shè)計參數(shù)樣本組合下懸臂梁加速度幅值。結(jié)合蒙特卡洛及隨機(jī)抽樣方法，選8萬個彈性模量E與材料密度ρ的參數(shù)樣本組合，計算懸臂梁加速度幅值。加速度響應(yīng)最大值非參數(shù)核密度估計［26］見圖5。由表2看出，在不同置信度條件下，加速度響應(yīng)最大值的置信區(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é)論

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

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

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

［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］吳立人.仿真模型有效性定量確認(rèn)方法［J］.導(dǎo)彈與航天運(yùn)載技術(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)動力學(xué)有限元模型確認(rèn)方法研究［J］.應(yīng)用力學(xué)學(xué)報，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］鄧小剛，宗文剛，張來平，等.計算流體力學(xué)中的驗證與確認(rèn)［J］.力學(xué)進(jìn)展，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］張保強(qiáng)，郭勤濤，陳國平.模型確認(rèn)熱傳導(dǎo)挑戰(zhàn)問題求解的貝葉斯方法［J］.航空學(xué)報，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īng)面方法在氣動優(yōu)化設(shè)計中的應(yīng)用［J］.工程熱物理學(xué)報，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］許瑞飛，宋文萍，韓忠華.改進(jìn)Kriging模型在翼型氣動優(yōu)化設(shè)計中的應(yīng)用研究［J］.西北工業(yè)大學(xué)學(xué)報，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］.機(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在邊坡可靠度分析中的應(yīng)用［J］.巖土力學(xué)，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ì)濤，劉 飛，張為華，等.基于拉丁超立方抽樣及響應(yīng)面的結(jié)構(gòu)模糊分析［J］.機(jī)械強(qiáng)度，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)計分析與應(yīng)用:40個案例分析［M］.北京:北京航空航天大學(xué)出版社，2010:355-369.