基于Kriging代理模型的中頻結(jié)構(gòu)振動(dòng)可靠性分析

方琛 秦朝紅 任方 張亞輝

(1 大連理工大學(xué) 工程力學(xué)系 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室 大連 116024;2 北京強(qiáng)度環(huán)境研究所 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室 北京 100076)

0 引言

工程結(jié)構(gòu)在制造、運(yùn)輸、裝配及服役過程不可避免地存在各種不確定性因素，所以在動(dòng)態(tài)激勵(lì)作用下單純地計(jì)算結(jié)構(gòu)的動(dòng)力學(xué)響應(yīng)或者響應(yīng)的極值會(huì)忽視不確定性的影響。隨著計(jì)算能力和精細(xì)化建模需求的提升，結(jié)合概率論準(zhǔn)確描述不確定性因素對(duì)動(dòng)力學(xué)響應(yīng)的影響，開展結(jié)構(gòu)的動(dòng)力學(xué)可靠度分析受到了廣泛關(guān)注。低頻問題在采用諸如有限元（Finite Element Method，F(xiàn)EM）等數(shù)值方法實(shí)現(xiàn)模型化后，通過概率模型結(jié)合抽樣或矩方法等即可得到結(jié)構(gòu)的失效概率[1]。隨著頻率的升高，繼續(xù)采用FEM建模需要將網(wǎng)格劃分得足夠細(xì)，這會(huì)導(dǎo)致龐大的自由度數(shù)。在計(jì)算量難以負(fù)擔(dān)的同時(shí)，還會(huì)出現(xiàn)參數(shù)敏感性問題。當(dāng)頻率足夠高時(shí)，充分考慮不確定性的集總參數(shù)的統(tǒng)計(jì)能量分析[2]（Statistical Energy Analysis，SEA）應(yīng)運(yùn)而生?？梢酝ㄟ^響應(yīng)的方差及近似概率分布描述非參數(shù)隨機(jī)的SEA模型的離散性[2-3]。

然而，實(shí)際工程結(jié)構(gòu)在寬頻激勵(lì)下，不同組件的振動(dòng)波長、模態(tài)密度等動(dòng)力學(xué)特性差異較大，就會(huì)表現(xiàn)出中頻振動(dòng)行為。傳統(tǒng)的低頻或高頻建模方法都無法直接應(yīng)用。目前主要有兩種解決方案，一是傳統(tǒng)方法向中頻段的擴(kuò)展[4-5]，二是基于子系統(tǒng)思想分別采用不同建模方法得到的中頻混合模型。其中Shorter和Langley基于擴(kuò)散場互易關(guān)系[6]提出的混合有限元-統(tǒng)計(jì)能量分析方法（Hybrid Finite Element-Statistical Energy Analysis method，簡稱混合FE-SEA方法）[7]理論最為完善，在諸如航空航天等領(lǐng)域應(yīng)用廣泛[8-9]。后續(xù)學(xué)者嘗試采用不同的方法替換FEM和SEA建立了不同的混合模型[10],[11]。為了實(shí)現(xiàn)中頻結(jié)構(gòu)振動(dòng)的可靠性分析，本文將采用混合FE-SEA方法實(shí)現(xiàn)模型化。

需要強(qiáng)調(diào)的是，混合FE-SEA方法將短波組件定義為統(tǒng)計(jì)性子系統(tǒng)，采用了集總參數(shù)的SEA建模，計(jì)算得到響應(yīng)是已經(jīng)考慮了短波組件非參數(shù)不確定性的。而長波組件是完全確定性的，與實(shí)際不符。隨著現(xiàn)在計(jì)算能力的發(fā)展和對(duì)不確定性研究的需求，安全系數(shù)法已經(jīng)不能滿足現(xiàn)有結(jié)構(gòu)的設(shè)計(jì)分析需求。Cicirello基于混合FE-SEA方法首次開展了中頻振動(dòng)的結(jié)構(gòu)可靠性分析[12]。Cicirello分別采用蒙特卡洛模擬（Monte-Carlo Simulation，MCS）和漸進(jìn)近似拉普拉斯方法計(jì)算得到了彈簧-板組合系統(tǒng)的中頻結(jié)構(gòu)可靠度。

混合FE-SEA方法的計(jì)算受制于復(fù)數(shù)總動(dòng)剛度矩陣的求逆運(yùn)算及譜分析形式，同等自由度規(guī)模下計(jì)算負(fù)擔(dān)遠(yuǎn)大于FEM。那么對(duì)混合模型直接采用MCS求解可靠度的將進(jìn)一步加劇計(jì)算負(fù)擔(dān)。此外基于混合FE-SEA方法建立的中頻問題的功能函數(shù)也格外復(fù)雜，采用傳統(tǒng)的可靠性分析方法求解比較困難。本文采用Kriging代理模型[13]-[15]直接建立隨機(jī)參數(shù)與響應(yīng)的關(guān)系，從而避免復(fù)雜的原物理模型重分析。Kriging模型的精度依賴于樣本點(diǎn)的選取，對(duì)于混合方法建立的復(fù)雜的中頻振動(dòng)求解模型，需要選取大量的樣本點(diǎn)，這樣會(huì)產(chǎn)生較大的計(jì)算量。為了在減少樣本數(shù)量的前提下保證計(jì)算精度，可以采用拉丁超立方試驗(yàn)設(shè)計(jì)[16]（Latin Hypercube Sampling，LHS）得到布滿樣本空間的試驗(yàn)組合，實(shí)現(xiàn)Kriging代理模型的全滿空間設(shè)計(jì)。故本文主要工作是基于混合FE-SEA方法的模型化，進(jìn)一步考慮隨機(jī)參數(shù)構(gòu)建了可靠性分析的功能函數(shù)，采用LHS結(jié)合Kriging模型構(gòu)建代理模型實(shí)現(xiàn)中頻振動(dòng)的結(jié)構(gòu)可靠度分析。

1 中頻振動(dòng)的混合FE-SEA模型

要實(shí)現(xiàn)中頻振動(dòng)的結(jié)構(gòu)可靠性分析，首先需要建立中頻振動(dòng)的響應(yīng)分析模型。本文采用經(jīng)典的混合FE-SEA方法實(shí)現(xiàn)中頻振動(dòng)的模型化?；旌螰E-SEA方法將復(fù)雜系統(tǒng)的各個(gè)組件區(qū)分為長波變形子系統(tǒng)和短波變形子系統(tǒng)，然后分別采用FEM和SEA建模，最后通過擴(kuò)散場互易關(guān)系實(shí)現(xiàn)系統(tǒng)響應(yīng)的耦合求解。本節(jié)將對(duì)中頻結(jié)構(gòu)的混合FE-SEA方法的基本理論進(jìn)行介紹。

1.1 總運(yùn)動(dòng)方程

如圖1所示，混合FE-SEA方法將短波組件定義為統(tǒng)計(jì)性子系統(tǒng)，認(rèn)為其是“直接場”與“混響場”的疊加。然后將混響場對(duì)直接場確定性邊界的作用定義為“混響場荷載”[6]。

圖1 混合FE-SEA模型的耦合示意圖 Fig. 1 Coupling of the hybrid FE-SEA method

為開展可靠性分析，本文進(jìn)一步通過概率模型考慮長波組件的參數(shù)隨機(jī)特性，并將長波組件定義為隨機(jī)子系統(tǒng)。首先建立系統(tǒng)的總運(yùn)動(dòng)方程，其描述的是直接場效應(yīng)和長波組件，寫作

其中()H表示共軛轉(zhuǎn)置，E[]表示集合平均，nj為第j個(gè)統(tǒng)計(jì)性子系統(tǒng)的模態(tài)密度，ω為激勵(lì)頻率。

為了聯(lián)立式(2)所示的擴(kuò)散場互易關(guān)系，求總運(yùn)動(dòng)方程(1)建立的自由度q的互譜，表達(dá)式為

1.2 功率流平衡方程

根據(jù)功率流平衡關(guān)系，第j個(gè)統(tǒng)計(jì)性子系統(tǒng)的功率流平衡方程表示為[7]

hij表示第j個(gè)統(tǒng)計(jì)性子系統(tǒng)混響場的單位模態(tài)密度輸入到第i個(gè)統(tǒng)計(jì)性子系統(tǒng)直接場的功率，表達(dá)式為

綜合式(5)～(7)，將m個(gè)統(tǒng)計(jì)性子系統(tǒng)的功率流平衡方程聯(lián)立為矩陣形式，表示為

E是由m個(gè)統(tǒng)計(jì)性子系統(tǒng)的能量組成的未知向量，是隨機(jī)子系統(tǒng)對(duì)各統(tǒng)計(jì)性子系統(tǒng)的輸入功率，N是各統(tǒng)計(jì)性子系統(tǒng)的模態(tài)密度組成的對(duì)角矩陣。L為系數(shù)矩陣，該系數(shù)矩陣的表達(dá)式為

求解式(8)即可得到所有統(tǒng)計(jì)性子系統(tǒng)的能量，回代到式(3)即可求解得到的隨機(jī)子系統(tǒng)的響應(yīng)Sqq。從系數(shù)矩陣L的表達(dá)式可以發(fā)現(xiàn)，混合模型求解需要首先計(jì)算總動(dòng)剛度矩陣的逆，然后進(jìn)行大量的求和運(yùn)算，計(jì)算負(fù)擔(dān)很大。

2 中頻振動(dòng)的結(jié)構(gòu)可靠性分析

上一節(jié)通過混合FE-SEA方法實(shí)現(xiàn)了中頻振動(dòng)的模型化，隨后即可開展結(jié)構(gòu)可靠性分析。本節(jié)將建立可靠性分析的功能函數(shù)，并通過構(gòu)造Kriging代理模型以避免復(fù)雜且計(jì)算量龐大的原物理模型的求解。

2.1 功能函數(shù)列式

基于混合FE-SEA方法進(jìn)行結(jié)構(gòu)中頻振動(dòng)的模型化時(shí)，統(tǒng)計(jì)性子系統(tǒng)采用集總參數(shù)的SEA方法建模時(shí)充分考慮了非參數(shù)隨機(jī)性。在忽略了短波組件離散性基礎(chǔ)上，本文假設(shè)可靠性分析的n個(gè)隨機(jī)參數(shù)x=[x1,…xn]來源于定義為隨機(jī)子系統(tǒng)的長波組件的相關(guān)參數(shù)。將通過混合FE-SEA中頻隨機(jī)模型求解的響應(yīng)表示為Q(x)。進(jìn)一步可以將結(jié)構(gòu)的極限狀態(tài)通過功能函數(shù)來表示

G(x)由結(jié)構(gòu)的集合平均響應(yīng)Q(x)和響應(yīng)界限值B來決定。Z＞0、Z＜0和Z=0就表示結(jié)構(gòu)處于安全、失效和極限狀態(tài)。通過采樣方法計(jì)算可靠度表達(dá)式為

即可靠度為安全狀態(tài)的樣本數(shù)np(Z＞0)占總樣本數(shù)nsum的百分比。

2.2 Kriging代理模型

中頻問題的功能函數(shù)(10)由于響應(yīng)Q(x)與隨機(jī)參數(shù)x=[x1,…xn]之間運(yùn)算關(guān)系格外復(fù)雜，難以采用傳統(tǒng)的動(dòng)力學(xué)可靠度方法求解。同時(shí)采用混合FE-SEA方法建立1.2節(jié)的功率流平衡方程時(shí)，必須直接對(duì)總動(dòng)剛度進(jìn)行求逆，然后進(jìn)行一系列求和運(yùn)算。如果隨機(jī)子系統(tǒng)的自由度規(guī)模較大，則計(jì)算負(fù)擔(dān)將難以承受。故采用Kriging方法建立等效響應(yīng)求解的代理模型，結(jié)合滿樣本空間的LHS控制樣本點(diǎn)的選取，就可以建立式(10)所需的響應(yīng)求解的等效方程。

Kriging代理模型是一個(gè)基于線性回歸分析建立的改進(jìn)擬合模型。該模型包括回歸項(xiàng)和隨機(jī)項(xiàng)兩部分，對(duì)于一個(gè)函數(shù)y(x0)，其可以表示為

x0表示待測(cè)點(diǎn)；fT(x0)表示給定的多項(xiàng)式回歸模型；β表示待定的回歸系數(shù)；z(x0)表示隨機(jī)過程。z(x0)的均值為0，標(biāo)準(zhǔn)差為σz，協(xié)方差表示為

R(x1,x2)表示樣本點(diǎn)中任意兩點(diǎn)x1、x2之間與距離dj有關(guān)的變異函數(shù)。R(x1,x2)由核函數(shù)計(jì)算，有諸多相關(guān)方程[17]，廣泛使用的是高斯相關(guān)方程，寫作

θj代表各向異性參數(shù)；p∈[0,2]表示光滑程度參數(shù)。

通過樣本數(shù)據(jù)X和對(duì)應(yīng)響應(yīng)Y，可以得到相關(guān)函數(shù)矩陣R?？紤]待測(cè)點(diǎn)x0與N個(gè)樣本點(diǎn)X =[x1,…xN]的相關(guān)向量r (x0)，采用最大似然估計(jì)得到θj優(yōu)化值。進(jìn)一步采用最小二乘法可得最大似然估計(jì)因子

z(x)的方差估計(jì)表達(dá)式為

則Kriging代理模型在待測(cè)點(diǎn)x0處的預(yù)測(cè)值為

至此建立得到Kriging代理模型。

建立的Kriging代理模型的計(jì)算精度與效率依賴于樣本點(diǎn)的選取，為了有效的選擇充滿空間的少量樣本點(diǎn)，本文選擇比較主流的LHS[16]采樣方法。通過LHS建立滿空間的nk個(gè)樣本點(diǎn)，然后通過物理模型計(jì)算得到相應(yīng)響應(yīng)Q(x)。最后結(jié)合式(17)就可以建立Kriging代理模型，即一個(gè)近似的響應(yīng)分析模型。

首先，建立代理模型所需的LHS樣本點(diǎn)要遠(yuǎn)少于直接針對(duì)物理模型進(jìn)行MCS的樣本點(diǎn)數(shù)。其次代理模型的計(jì)算效率要遠(yuǎn)高于原物理模型，所以針對(duì)代理模型進(jìn)行MCS即可迅速得到響應(yīng)可靠度。本文采用成熟的DACE工具箱[18]建立Kriging模型，該工具箱簡便實(shí)用，被工程師和研究者們廣泛采用。

3 數(shù)值算例

本節(jié)將對(duì)如圖2所示的一個(gè)框架梁栓接板組合系統(tǒng)開展中頻結(jié)構(gòu)動(dòng)力學(xué)可靠度分析。該結(jié)構(gòu)由一個(gè)邊長0.7m為正方體框架梁和四個(gè)長1.1m、寬0.6m的矩形板組成。梁截面為矩形，各板通過四個(gè)點(diǎn)栓接在框架梁上，各栓接點(diǎn)距離框架梁頂點(diǎn)0.1m。在原點(diǎn)o作用一個(gè)沿y軸正向的200Hz的單頻點(diǎn)荷載?？蚣芰旱牟牧蠟殇?，阻尼損耗因子為0.01，板為鋁，阻尼損耗因子為0.02。

圖2 框架梁栓接板結(jié)構(gòu) Fig. 2 A framework-panel system

該模型其余材料屬性、幾何參數(shù)、荷載幅值等均為隨機(jī)參數(shù)，并考慮環(huán)境影響的不確定性。相關(guān)隨機(jī)參數(shù)的具體統(tǒng)計(jì)信息及概率分布見表1。

表1 隨機(jī)變量的概率分布參數(shù) Table 1 Probability distribution parameters of random variables

為開展可靠性分析，首先需要基于混合FE-SEA方法實(shí)現(xiàn)模型化。而混合模型的建模需要首先根據(jù)振動(dòng)波長進(jìn)行子系統(tǒng)劃分，200Hz激勵(lì)下梁和板的波長見表2，計(jì)算采用名義參數(shù)。如圖3 a)，將長波變形的梁和板面內(nèi)振動(dòng)定義為隨機(jī)子系統(tǒng)，采用FEM建模，網(wǎng)格尺寸為0.05m。如圖3 b)，采用SEA建模，將板短波變形的彎曲振動(dòng)定義為統(tǒng)計(jì)性子系統(tǒng)。

表2 振動(dòng)波長 Table 2 Vibrational wavelength

圖3 混合FE-SEA模型的劃分 Fig. 3 Subsystems of hybrid FE-SEA model

通過混合FE-SEA方法實(shí)現(xiàn)模型化后，即可開展可靠性分析。對(duì)于隨機(jī)子系統(tǒng)，選擇梁上的測(cè)點(diǎn)A(0.5, 0.7, 0.7)的速度幅值互譜Svv=1.66 ×10-5m2/s2為臨界變形。對(duì)于統(tǒng)計(jì)性子系統(tǒng)，選擇z=0.7m平面上的板的能量E4=5.42 ×10-4J 作為臨界能量。

對(duì)混合FE-SEA方法建立的原物理模型進(jìn)行nm=105個(gè)樣本點(diǎn)的MCS[12]，將其計(jì)算結(jié)果作為可靠度參考解。對(duì)物理模型選取數(shù)量nk=500,600，1000的LHS樣本（隨機(jī)參數(shù)置信區(qū)間不小于99%）建立Kriging代理模型，可靠度計(jì)算結(jié)果如表3所示。通過nk=500個(gè)樣本點(diǎn)建立的代理模型計(jì)算的可靠度已經(jīng)接近MCS的結(jié)果，誤差均小于2%。隨著LHS樣本點(diǎn)的增多，誤差逐漸減小到低于0.5%。值得注意的是，該混合模型單次采樣計(jì)算都非常耗時(shí)，其中大部分時(shí)間耗費(fèi)在總動(dòng)剛度矩陣的求逆及相關(guān)的求和運(yùn)算中。而采用代理模型計(jì)算近似響應(yīng)的時(shí)間可以忽略不計(jì)。所以采用LHS采樣建立Kriging代理模型，求解可靠度的方法通過較小的計(jì)算量實(shí)現(xiàn)了中頻振動(dòng)的結(jié)構(gòu)可靠度求解。即使是LHS采樣數(shù)nk=1000，代理模型的計(jì)算時(shí)間只是nm=105次MCS的1%而已。而LHS采樣數(shù)為nk=500,600次時(shí)間，計(jì)算效率可以在MCS基礎(chǔ)上分別提升200、167倍。

表3 結(jié)構(gòu)可靠度分析結(jié)果(%) Table 3 Reliability analysis results (%)

進(jìn)一步驗(yàn)證更高可靠度下的結(jié)果。將測(cè)點(diǎn)A(0.5, 0.7, 0.7)的速度幅值互譜和z=0.7m平面上的板能量對(duì)應(yīng)的臨界值重新給定為Svv=1.94 ×10-5m2/s2和E4=5.71 ×10-4J 。計(jì)算結(jié)果如表4所示，此時(shí)MCS得到的可靠度均為99%。采用nk= 1000個(gè)LHS樣本建立了Kriging代理模型，計(jì)算得到的可靠度誤差均小于1%。

表4 結(jié)構(gòu)可靠度分析結(jié)果(%) Table 4 Reliability analysis results(%)

4 結(jié)論

本文開展了中頻振動(dòng)的結(jié)構(gòu)可靠性分析?；旌螰E-SEA方法將短波組件采用SEA建模，本文在忽略短波組件的離散性后通過概率模型描述長波變形組件的隨機(jī)性來開展可靠性分析。通過滿空間的LHS選取樣本點(diǎn)即可建立Kriging代理模型。針對(duì)一個(gè)框架梁栓接板組合結(jié)構(gòu)計(jì)算可靠度。發(fā)現(xiàn)較之對(duì)物理模型通過大量采樣的直接MCS的方法，通過少量的LHS樣本點(diǎn)建立的代理模型可以迅速計(jì)算得到可靠度，計(jì)算效率很高。隨著LHS樣本點(diǎn)增多，誤差會(huì)逐漸減小。