本發(fā)明屬于模型確認(rèn)技術(shù)領(lǐng)域，尤其設(shè)計一種基于kriging響應(yīng)面的有限元模型確認(rèn)方法。

背景技術(shù)：

在現(xiàn)代工程結(jié)構(gòu)的研發(fā)和設(shè)計過程中以及重大裝備或工程設(shè)施在服役階段的安全性和可靠性評估方面，計算仿真正發(fā)揮著越來越重要的作用，特別是產(chǎn)品設(shè)計階段以及復(fù)雜系統(tǒng)很難進(jìn)行全系統(tǒng)試驗的情形下，就必須依靠計算仿真進(jìn)行評估，因此仿真模型的精度和置信度至關(guān)重要。模型確認(rèn)在工程界的最初概念由美國能源部提出，主要應(yīng)用于戰(zhàn)略武器存儲管理的可靠性評估和決策中，并于1998年由aiaa起草了計算流體動力學(xué)仿真模型驗證和確認(rèn)的指南，oberkampf對此進(jìn)行了系統(tǒng)總結(jié)，并綜述了機械工程領(lǐng)域仿真模型確認(rèn)的發(fā)展。jung將模型確認(rèn)的分層體系應(yīng)用到一般工程產(chǎn)品的開發(fā)中。在計算仿真過程中，由于各種不確定性的存在使得計算仿真的預(yù)測結(jié)果與試驗結(jié)果之間往往存在很大差異，因此模型修正技術(shù)也得到了廣泛的應(yīng)用，重點對不確定性對模型預(yù)測的影響和修正方法進(jìn)行討論。模型確認(rèn)在國內(nèi)也逐漸受到重視，模型確認(rèn)的概念最早由張令彌引入到國內(nèi)，郭勤濤等將模型確認(rèn)應(yīng)用到具體研究中，王瑞利和鄧小剛分別對模型確認(rèn)在計算機程序和流體動力學(xué)中的應(yīng)用進(jìn)行討論。劉信恩對模型確認(rèn)中的貝葉斯框架進(jìn)行討論和簡化。然而，總體來看，模型確認(rèn)的基本流程尚不明確，還處在研究階段。

對于復(fù)雜結(jié)構(gòu)的輸入?yún)?shù)和輸出相應(yīng)特征之間的轉(zhuǎn)換關(guān)系，可以通過構(gòu)建代理模型來代替有限元模型，可以在短時間進(jìn)行大量隨機抽樣，大大減少計算量。近年來，響應(yīng)面模型、kriging模型和神經(jīng)網(wǎng)絡(luò)等代理模型快速發(fā)展。然而各種代理模型有其各自的優(yōu)缺點。響應(yīng)面模型易于實現(xiàn)，但是逼近非線性問題能力較差；kriging模型對非線性問題有較高的準(zhǔn)確度，但模型的獲取和使用難度較大，神經(jīng)網(wǎng)絡(luò)模型的精度和泛化能力依賴其網(wǎng)絡(luò)結(jié)構(gòu)和大量的學(xué)習(xí)樣本，有時存在“過學(xué)習(xí)”現(xiàn)象。同時，由于復(fù)雜模型的不確定性參數(shù)眾多，如何選取對目標(biāo)響應(yīng)相關(guān)性最大的參數(shù)也是研究的熱點。

技術(shù)實現(xiàn)要素：

本發(fā)明的發(fā)明目的是：為了解決以上問題，本發(fā)明提出了一種基于kriging響應(yīng)面的有限元模型確認(rèn)方法，為復(fù)雜系統(tǒng)計算仿真中的模型確認(rèn)工作提供便利，同時為產(chǎn)品的結(jié)構(gòu)設(shè)計和健康評估等提供依據(jù)。

本發(fā)明的技術(shù)方案是：一種基于kriging響應(yīng)面的有限元模型確認(rèn)方法，包括以下步驟：

a、對壓力容器模型進(jìn)行子系統(tǒng)劃分，建立有限元模型；

b、對步驟a中有限元模型的模型參數(shù)進(jìn)行參數(shù)相關(guān)性和參數(shù)靈敏度分析；

c、選取與最大等效應(yīng)力相關(guān)性最大的參數(shù)，采用拉丁超立方實驗設(shè)計方法構(gòu)建壓力容器模型最大等效應(yīng)力的kriging響應(yīng)面；

d、對步驟c中kriging響應(yīng)面進(jìn)行隨機抽樣，采用核密度估計方法對抽樣樣本概率分布和置信區(qū)間進(jìn)行分析；

e、將單次有限元計算結(jié)果與有限元仿真結(jié)果進(jìn)行對比，實現(xiàn)有限元模型確認(rèn)。

進(jìn)一步地，所述步驟b對有限元模型的模型參數(shù)進(jìn)行參數(shù)相關(guān)性和參數(shù)靈敏度分析具體為：

將有限元模型的模型參數(shù)分為控制參數(shù)和幾何材料參數(shù)，采用實驗設(shè)計模塊計算幾何材料參數(shù)對最大等效應(yīng)力和最大變形量的靈敏度，得到模型參數(shù)靈敏度柱狀圖；采用實驗設(shè)計模塊計算幾何材料參數(shù)對最大等效應(yīng)力的相關(guān)性，得到模型參數(shù)相關(guān)性矩陣。

進(jìn)一步地，所述步驟c中壓力容器模型最大等效應(yīng)力的kriging響應(yīng)面的數(shù)學(xué)模型具體為：

其中，為kriging響應(yīng)面在樣本點x處的響應(yīng)預(yù)測值，為估計后的基函數(shù)系數(shù)，f(x)為擬合函數(shù)的基函數(shù)，r為相關(guān)系數(shù)矩陣，yd為樣本的觀測值，f為f(x)的觀測值。

進(jìn)一步地，所述步驟d對kriging響應(yīng)面進(jìn)行隨機抽樣，采用核密度估計方法對抽樣樣本概率分布和置信區(qū)間進(jìn)行分析具體為：

采用蒙特卡洛及隨機抽樣方法，選取5000個半徑r、厚度th的參數(shù)樣本組合，對kriging響應(yīng)面進(jìn)行抽樣得到5000組最大等效應(yīng)力幅值，利用直方圖法對最大等效應(yīng)力進(jìn)行分析；采用非參數(shù)核密度估計方法計算等效應(yīng)力最大值的概率密度函數(shù)，通過設(shè)置不同的置信度水平，得到壓力容器模型的最大等效應(yīng)力的置信區(qū)間。

進(jìn)一步地，所述步驟e將單次有限元計算結(jié)果與有限元仿真結(jié)果進(jìn)行對比，實現(xiàn)有限元模型確認(rèn)具體為：

根據(jù)壓力容器模型設(shè)計低水平和高水平大氣壓力試驗，構(gòu)建兩組kriging響應(yīng)面，對每個kriging響應(yīng)面進(jìn)行抽樣得到認(rèn)證樣本，將兩組得到的認(rèn)證樣本利用核密度估計方法進(jìn)行擬合，得到kde曲線；將kde曲線與兩種實驗條件下的有限元仿真得到的最大等效應(yīng)力值進(jìn)行比較，實現(xiàn)有限元模型確認(rèn)。

本發(fā)明的有益效果是：本發(fā)明通過對壓力容器模型系統(tǒng)進(jìn)行分層設(shè)計，建立參數(shù)化有限元模型，并對模型參數(shù)進(jìn)行參數(shù)相關(guān)性和參數(shù)靈敏度分析，選取與目標(biāo)響應(yīng)相關(guān)性最大的參數(shù)構(gòu)建kriging響應(yīng)面，對響應(yīng)面隨機抽樣實現(xiàn)對目標(biāo)響應(yīng)的預(yù)測，利用核密度估計方法對抽樣樣本概率分布和置信區(qū)間進(jìn)行分析，并與單次有限元計算結(jié)果進(jìn)行對比實現(xiàn)模型確認(rèn)，從而為模型確認(rèn)提供便利，同時為產(chǎn)品的結(jié)構(gòu)設(shè)計和健康評估等提供依據(jù)。

附圖說明

圖1是本發(fā)明的基于kriging響應(yīng)面的有限元模型確認(rèn)方法的流程示意圖。

圖2是本發(fā)明實施例中壓力容器模型確認(rèn)分層示意圖。

圖3是本發(fā)明實施例中壓力容器模型試驗設(shè)計載荷施加示意圖。

圖4是本發(fā)明實施例中參數(shù)靈敏度柱狀圖。

圖5是本發(fā)明實施例中參數(shù)相關(guān)性矩陣圖。

圖6是本發(fā)明實施例中壓力容器模型最大等效應(yīng)力kriging響應(yīng)面示意圖。

圖7是本發(fā)明實施例中壓力容器模型最大等效應(yīng)力核密度估計示意圖。

圖8是本發(fā)明實施例中認(rèn)證實驗一的kriging響應(yīng)面示意圖。

圖9是本發(fā)明實施例中認(rèn)證實驗一的kde曲線示意圖。

圖10是本發(fā)明實施例中認(rèn)證實驗一的cdf曲線示意圖。

圖11是本發(fā)明實施例中認(rèn)證實驗二的kriging響應(yīng)面示意圖。

圖12是本發(fā)明實施例中認(rèn)證實驗二的kde曲線示意圖。

圖13是本發(fā)明實施例中認(rèn)證實驗二的cdf曲線示意圖。

具體實施方式

為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點更加清楚明白，以下結(jié)合附圖及實施例，對本發(fā)明進(jìn)行進(jìn)一步詳細(xì)說明。應(yīng)當(dāng)理解，此處所描述的具體實施例僅用以解釋本發(fā)明，并不用于限定本發(fā)明。

如圖1所示，為本發(fā)明的基于kriging響應(yīng)面的有限元模型確認(rèn)方法的流程示意圖。一種基于kriging響應(yīng)面的有限元模型確認(rèn)方法，包括以下步驟：

a、對壓力容器模型進(jìn)行子系統(tǒng)劃分，建立有限元模型；

b、對步驟a中有限元模型的模型參數(shù)進(jìn)行參數(shù)相關(guān)性和參數(shù)靈敏度分析；

c、選取與最大等效應(yīng)力相關(guān)性最大的參數(shù)，采用拉丁超立方實驗設(shè)計方法構(gòu)建壓力容器模型最大等效應(yīng)力的kriging響應(yīng)面；

d、對步驟c中kriging響應(yīng)面進(jìn)行隨機抽樣，采用核密度估計方法對抽樣樣本概率分布和置信區(qū)間進(jìn)行分析；

e、將單次有限元計算結(jié)果與有限元仿真結(jié)果進(jìn)行對比，實現(xiàn)有限元模型確認(rèn)。

在步驟a中，由于壓力容器模型存在固液氣耦合問題，不確定參數(shù)較多，因此本發(fā)明對壓力容器模型進(jìn)行子系統(tǒng)劃分，并建立相關(guān)的有限元模型；劃分后的子系統(tǒng)擁有較少的不確定參數(shù)，在分析過程中可以大大減少不確定參數(shù)的誤差耦合產(chǎn)生的影響。本發(fā)明通過對壓力容器的實際載荷和約束情況進(jìn)行分析，從環(huán)境復(fù)雜度和系統(tǒng)復(fù)雜度兩個方面進(jìn)行考慮，對仿真模型進(jìn)行子系統(tǒng)劃分。如圖2所示，為本發(fā)明實施例中壓力容器模型確認(rèn)分層示意圖。

本發(fā)明中壓力容器為一個密閉容器，由一個圓柱形殼體和兩個半球形殼體構(gòu)成，在圓柱面和球面連接端施加固定約束，壓力容器內(nèi)有一定的大氣壓強和一定高度的特殊液體。在pro/e軟件中建立壓力容器的參數(shù)化三維實體模型，并通過與ansysworkbench14.0的無縫鏈接建立有限元分析模型。

本發(fā)明根據(jù)有限元模型確認(rèn)流程和分層思想，對壓力容器模型進(jìn)行子系統(tǒng)劃分和試驗設(shè)計。試驗設(shè)計分為校準(zhǔn)試驗、確認(rèn)試驗和認(rèn)證試驗，分別施加不同層次的載荷進(jìn)行試驗。如圖3所示，為本發(fā)明實施例中壓力容器模型試驗設(shè)計載荷施加示意圖，其中左為校準(zhǔn)試驗，中為確認(rèn)試驗，右為認(rèn)證試驗。

(1)校準(zhǔn)試驗：常壓；

(2)確認(rèn)試驗：施加一定壓強p載荷；

(3)認(rèn)證試驗：施加一定壓強p和一定高度h液體混合載荷。

校準(zhǔn)試驗是通過試驗對子結(jié)構(gòu)模型參數(shù)(如材料彈性模量、泊松比、密度、比重等)進(jìn)行直接或間接測量，對參數(shù)的誤差進(jìn)行估計；確認(rèn)試驗通過構(gòu)建kriging響應(yīng)面模型，并對目標(biāo)響應(yīng)的概率分布和置信區(qū)間進(jìn)行分析；認(rèn)證試驗在確認(rèn)過的子系統(tǒng)模型基礎(chǔ)上，構(gòu)建認(rèn)證試驗有限元模型，并將有限元模型仿真結(jié)果與響應(yīng)面預(yù)測結(jié)果進(jìn)行比較。

在步驟b中，本發(fā)明將有限元模型的模型參數(shù)分為控制參數(shù)和幾何材料參數(shù)?？刂茀?shù)包括壓強、液體高度、液體比重，幾何材料參數(shù)包括長、半徑、厚度、彈性模量等。

利用ansysworkbench中的實驗設(shè)計模塊(doe)對確認(rèn)試驗中的不確定參數(shù)l、r、th、e進(jìn)行靈敏度分析，將參數(shù)l、r、th、e作為輸入?yún)?shù)，將最大等效應(yīng)力和最大變形量作為輸出參數(shù)。假設(shè)設(shè)結(jié)構(gòu)設(shè)計參數(shù)為x～n(μ,σ)的正態(tài)分布，μ為各設(shè)計參數(shù)初值，σ＝μα為變量標(biāo)準(zhǔn)差，α為設(shè)計參數(shù)的變異系數(shù)，五個設(shè)計參數(shù)的變異系數(shù)均取5％。如表1所示為壓力容器模型進(jìn)行doe分析的不確定性參數(shù)。

表1壓力容器模型不確定性參數(shù)取值

利用doe模塊產(chǎn)生25組設(shè)計點，并且對每個設(shè)計點進(jìn)行求解，得到25組壓力容器的最大等效應(yīng)力和最大總變形量的樣本點。如圖4所示，為本發(fā)明實施例中參數(shù)靈敏度柱狀圖，可以看出半徑r和厚度th對最大應(yīng)力的靈敏度性是最高的，模型的長度l和彈性模量e的靈敏度則相當(dāng)小，基本可以忽略不計；而材料參數(shù)彈性模量e對最大應(yīng)變量的靈敏度是最大的。

本發(fā)明利用workbench的doe模塊進(jìn)行實驗設(shè)計可以得到詳細(xì)的輸入?yún)?shù)與輸出參數(shù)的具體相關(guān)性數(shù)值。在參數(shù)相關(guān)性分析當(dāng)中，若輸入?yún)?shù)對輸出參數(shù)的相關(guān)度大于80％，則表示相關(guān)性較高，反之，若相關(guān)度小于10％，則表示相關(guān)性較低。如圖5所示為本發(fā)明實施例中參數(shù)相關(guān)性矩陣圖，如表2所示為參數(shù)相關(guān)性矩陣表，可以發(fā)現(xiàn)輸入?yún)?shù)中的壓力容器的厚度th與最大等效應(yīng)力的相關(guān)度達(dá)到了81.83％，半徑r與最大等效應(yīng)力的相關(guān)度為54.93％，其余輸入?yún)?shù)對最大等效應(yīng)力的相關(guān)度均小于10％，可以認(rèn)為壓力容器的半徑和厚度對最大等效應(yīng)力的相關(guān)性相對較高，壓力容器的長度和材料彈性模量對最大等效應(yīng)力的相關(guān)性很低。利用相關(guān)性分析可以選取與關(guān)心量(最大等效應(yīng)力)相關(guān)性最大的參數(shù)(半徑r和厚度th)來構(gòu)建響應(yīng)面。

表2參數(shù)相關(guān)性矩陣表

在步驟c中，本發(fā)明采用拉丁超立方試驗設(shè)計方法(latinhypercubedesign,ld)來構(gòu)建kriging響應(yīng)面。選取了25組設(shè)計樣本點，將輸入?yún)?shù)的x軸為壓力容器半徑r，輸入?yún)?shù)的y軸為壓力容器殼體厚度th，z軸為最大等效應(yīng)力(maxeqv.stress)。如圖6所示為本發(fā)明實施例中壓力容器模型最大等效應(yīng)力kriging響應(yīng)面示意圖

本發(fā)明的kriging響應(yīng)面模型以多項式逼近的方式對計算仿真模型進(jìn)行近似描述，其表達(dá)式為：

其中，fj(x)是擬合函數(shù)的基函數(shù)，βj是基函數(shù)的系數(shù)，z(x)是用于擬合的偏差函數(shù)。kriging插值方法一般認(rèn)為不同插值點處的擬合偏差量不是相互獨立的，并假設(shè)偏差函數(shù)是一種隨機過程z(x)，隨機過程的均值為0，方差為σ²，并且協(xié)方差不為0。隨機兩點t和u處的協(xié)方差函數(shù)定義為：

cov[z(t),z(u)]＝σ²ρ(t,u；θ)

其中，ρ(t,u；θ)為兩點間的相關(guān)函數(shù)；θ為相關(guān)函數(shù)的參數(shù)，該參數(shù)用來衡量兩個樣本點t與u之間的相關(guān)性隨著兩個樣本點間距增加的衰減度，此參數(shù)越小，則所構(gòu)建的響應(yīng)面就越光滑。

在確定相關(guān)函數(shù)ρ和樣本的觀測值yd＝[y(x1),y(x2),…y(xn)]^t后，還要依據(jù)樣本觀察值計算隨機過程方差σ²、基函數(shù)系數(shù)β以及相關(guān)函數(shù)參數(shù)θ，將估計后的三個參數(shù)記作即為kriging模型的超參數(shù)(hyperparameters)。利用最大似然估計法構(gòu)建觀測樣本值的聯(lián)合概率分布函數(shù)p(yd)，通過將聯(lián)合概率分布函數(shù)最大化來得到超參數(shù)具體數(shù)值：

l(β,σ²,θ)＝p(yd)

其中，

f＝[f(x1),f(x2),…f(xn)]^t

r為相關(guān)系數(shù)矩陣，rij＝ρ[xi,xi；θ],(1≤i,j≤n)。

參數(shù)β、θ表達(dá)式為

得到：

相關(guān)系數(shù)θ利用數(shù)值方法進(jìn)行求解，對極大似然函數(shù)求極值得到其數(shù)值的大小。記r(x)＝[ρ(x,x1),ρ(x,x2),…,ρ(x,xn),]^t，則kriging響應(yīng)面模型在任意設(shè)計樣本點x處的響應(yīng)預(yù)測表達(dá)式為：

為了對比不同的響應(yīng)面實驗設(shè)計方法對響應(yīng)面的預(yù)測精度的影響，在相同樣本點數(shù)(25組)的條件下，采用另外兩種實驗設(shè)計方法(水平全因子設(shè)計和復(fù)合中心設(shè)計)來構(gòu)建kriging響應(yīng)面模型。同時，利用響應(yīng)面精度檢驗指標(biāo)，對所構(gòu)建的三種響應(yīng)面精度進(jìn)行檢驗。本發(fā)明利用5個隨機樣本點作為檢驗點，如表3所示為kriging響應(yīng)面精度檢驗結(jié)果。

表3kriging響應(yīng)面精度檢驗結(jié)果

對比三種試驗設(shè)計方法構(gòu)建的kriging響應(yīng)面精度檢驗指標(biāo)r²、rmse可以發(fā)現(xiàn)，拉丁超立方設(shè)計方法可以在相同數(shù)量的樣本點條件下，構(gòu)建更加精確的kriging響應(yīng)面模型。

綜上分析，本發(fā)明利用拉丁超立方實驗設(shè)計方法來構(gòu)建壓力容器最大等效應(yīng)力的kriging響應(yīng)面是合理有效的。

在步驟d中，本發(fā)明的kriging響應(yīng)面以試驗設(shè)計參數(shù)樣本點為基礎(chǔ)，來對非試驗設(shè)計參數(shù)樣本組合下壓力容器殼體表面最大等效應(yīng)力幅值進(jìn)行預(yù)測。結(jié)合蒙特卡洛及隨機抽樣方法對響應(yīng)面進(jìn)行抽樣，選5000個半徑r、厚度th的參數(shù)樣本組合，對響應(yīng)面進(jìn)行抽樣得到5000組最大等效應(yīng)力幅值。首先利用直方圖法對最大等效應(yīng)力進(jìn)行分析，觀察其分布情況。

由于壓力容器最大等效應(yīng)力的分布類型未知，采用非參數(shù)核密度估計方法(kerneldensityestimation,kde)計算得到等效應(yīng)力最大值的概率密度函數(shù)，如圖7所示為本發(fā)明實施例中壓力容器模型最大等效應(yīng)力核密度估計示意圖；并通過設(shè)置不同的置信度水平，得到壓力容器的最大等效應(yīng)力的置信區(qū)間，如表4所示為等效應(yīng)力幅值在不同置信度下的置信區(qū)間，可以發(fā)現(xiàn)在不同置信度條件下，等效應(yīng)力最大值的上下限并未隨置信度的下降發(fā)生顯著變化。

表4等效應(yīng)力幅值在不同置信度下的置信區(qū)間

在步驟e中，本發(fā)明在建立的有限元模型中加入實驗特例對建立的有限元模型進(jìn)行認(rèn)證評估。本發(fā)明設(shè)計低水平和高水平大氣壓力的試驗來對確認(rèn)試驗進(jìn)行認(rèn)證評估，選取壓力容器測試數(shù)據(jù)中的兩組數(shù)據(jù)，分別設(shè)計認(rèn)證實驗一：p＝36.725psi，χ＝0.5，h＝13in；認(rèn)證實驗二：p＝52.884psi，χ＝0.7，h＝45in進(jìn)行試驗，依照確認(rèn)試驗的流程得到在兩組控制參數(shù)下的kriging響應(yīng)面，對響應(yīng)面進(jìn)行蒙特卡洛隨機抽樣后得到頻率累積直方圖、核密度估計曲線和累積分布函數(shù)曲線。

從確認(rèn)過的子結(jié)構(gòu)進(jìn)行實驗設(shè)計，在兩組試驗條件下，構(gòu)建兩組kriging響應(yīng)面，對每個響應(yīng)面代理模型進(jìn)行抽樣5000次，就得到5000個認(rèn)證樣本，將兩組得到的認(rèn)證樣本利用核密度估計方法進(jìn)行擬合，得到kde曲線。同時，將兩種實驗條件下的有限元仿真得到的最大等效應(yīng)力值與kde曲線進(jìn)行對比。如圖8所示為本發(fā)明實施例中認(rèn)證實驗一的kriging響應(yīng)面示意圖。如圖9所示為本發(fā)明實施例中認(rèn)證實驗一的kde曲線示意圖。如圖10所示為本發(fā)明實施例中認(rèn)證實驗一的cdf曲線示意圖。如圖11所示為本發(fā)明實施例中認(rèn)證實驗二的kriging響應(yīng)面示意圖。如圖12所示為本發(fā)明實施例中認(rèn)證實驗二的kde曲線示意圖。如圖13所示為本發(fā)明實施例中認(rèn)證實驗二的cdf曲線示意圖。通過kde曲線和有限元仿真得到的最大等效應(yīng)力值進(jìn)行比較可知，響應(yīng)面預(yù)測數(shù)據(jù)與仿真數(shù)據(jù)有一定的偏差，但相差不太大。大部分的響應(yīng)面預(yù)測結(jié)果都要比有限元仿真結(jié)果要小，因此利用kriging響應(yīng)面代理模型進(jìn)行預(yù)測的結(jié)果偏保守，從而實現(xiàn)有限元模型確認(rèn)。

本領(lǐng)域的普通技術(shù)人員將會意識到，這里所述的實施例是為了幫助讀者理解本發(fā)明的原理，應(yīng)被理解為本發(fā)明的保護范圍并不局限于這樣的特別陳述和實施例。本領(lǐng)域的普通技術(shù)人員可以根據(jù)本發(fā)明公開的這些技術(shù)啟示做出各種不脫離本發(fā)明實質(zhì)的其它各種具體變形和組合，這些變形和組合仍然在本發(fā)明的保護范圍內(nèi)。