方琛 秦朝紅 任方 張亞輝

(1 大連理工大學(xué) 工程力學(xué)系 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(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)波長(zhǎ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ò)散場(chǎng)互易關(guān)系[6]提出的混合有限元-統(tǒng)計(jì)能量分析方法（Hybrid Finite Element-Statistical Energy Analysis method，簡(jiǎn)稱混合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ù)不確定性的。而長(zhǎng)波組件是完全確定性的，與實(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)的模型化。混合FE-SEA方法將復(fù)雜系統(tǒng)的各個(gè)組件區(qū)分為長(zhǎng)波變形子系統(tǒng)和短波變形子系統(tǒng)，然后分別采用FEM和SEA建模，最后通過擴(kuò)散場(chǎng)互易關(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)為其是“直接場(chǎng)”與“混響場(chǎng)”的疊加。然后將混響場(chǎng)對(duì)直接場(chǎng)確定性邊界的作用定義為“混響場(chǎng)荷載”[6]。

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

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

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

為了聯(lián)立式(2)所示的擴(kuò)散場(chǎng)互易關(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)混響場(chǎng)的單位模態(tài)密度輸入到第i個(gè)統(tǒng)計(jì)性子系統(tǒng)直接場(chǎ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)的長(zhǎ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模型，該工具箱簡(jiǎn)便實(shí)用，被工程師和研究者們廣泛采用。

3 數(shù)值算例

本節(jié)將對(duì)如圖2所示的一個(gè)框架梁栓接板組合系統(tǒng)開展中頻結(jié)構(gòu)動(dòng)力學(xué)可靠度分析。該結(jié)構(gòu)由一個(gè)邊長(zhǎng)0.7m為正方體框架梁和四個(gè)長(zhǎng)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)波長(zhǎng)進(jìn)行子系統(tǒng)劃分，200Hz激勵(lì)下梁和板的波長(zhǎng)見表2，計(jì)算采用名義參數(shù)。如圖3 a)，將長(zhǎng)波變形的梁和板面內(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)波長(zhǎ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建模，本文在忽略短波組件的離散性后通過概率模型描述長(zhǎng)波變形組件的隨機(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ì)逐漸減小。