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.