基于響應(yīng)面與靈敏度分析的區(qū)間不確定性參數(shù)識(shí)別方法

陳學(xué)前， 沈展鵬， 劉信恩

(1. 中國(guó)工程物理研究院 總體工程研究所，四川 綿陽(yáng) 621999;2. 工程材料與結(jié)構(gòu)沖擊振動(dòng)四川省重點(diǎn)實(shí)驗(yàn)室，四川 綿陽(yáng) 621999)

可靠的動(dòng)力學(xué)模型對(duì)于結(jié)構(gòu)的動(dòng)力響應(yīng)預(yù)測(cè)和結(jié)構(gòu)動(dòng)態(tài)設(shè)計(jì)都具有十分重要的意義。因此，常常需要根據(jù)結(jié)構(gòu)振動(dòng)或模態(tài)試驗(yàn)結(jié)果對(duì)結(jié)構(gòu)的有限元模型進(jìn)行修正，使有限元分析結(jié)果能更好服務(wù)結(jié)構(gòu)的工程設(shè)計(jì)。為此，有限元建模與修正技術(shù)在過(guò)去30多年取得了長(zhǎng)足的發(fā)展[1-6]，但當(dāng)前絕大多數(shù)的模型修正方法都屬于確定性方法，沒(méi)有考慮結(jié)構(gòu)參數(shù)和響應(yīng)的不確定性，大大降低了其工程實(shí)際的應(yīng)用價(jià)值。然而，由于工程問(wèn)題普遍存在不確定性，如結(jié)構(gòu)材料參數(shù)與幾何尺寸的不確定性，結(jié)構(gòu)在服役期間的各類裝配帶來(lái)的裝配不確定性，試驗(yàn)測(cè)試數(shù)據(jù)的不確定性等，迫切需要開(kāi)展考慮不確定性的結(jié)構(gòu)動(dòng)力學(xué)有限元模型修正方法研究?？紤]不確定性的有限元模型修正可以分為概率模型修正和區(qū)間模型修正兩大類。在概率模型修正方面，最初的概率模型修正是采用統(tǒng)計(jì)方法僅考慮噪聲引起的試驗(yàn)結(jié)果隨機(jī)性的問(wèn)題[7]，也有研究者將貝葉斯方法應(yīng)用于模型修正以獲得更加可靠的識(shí)別結(jié)果[8-9]，近年來(lái)，為提高模型修正效率，有學(xué)者結(jié)合攝動(dòng)法與靈敏度分析方法開(kāi)展了不確定性模型修正[10-13]。但概率修正方法的精度依賴于對(duì)結(jié)構(gòu)參數(shù)和響應(yīng)概率分布特性的準(zhǔn)確估計(jì)，這類方法需要大量的試驗(yàn)數(shù)據(jù)以建立準(zhǔn)確的概率分布函數(shù)，增大了工程應(yīng)用的難度。相比較，采用區(qū)間分析的不確定性有限元模型修正對(duì)試驗(yàn)數(shù)據(jù)沒(méi)有嚴(yán)格要求，為此，當(dāng)前許多學(xué)者提出了區(qū)間響應(yīng)面模型以及模型修正的區(qū)間反演方法[14-16]。但是，由于區(qū)間數(shù)運(yùn)算法則和傳統(tǒng)數(shù)學(xué)算法非常不同，采用區(qū)間分析方法，待修正參數(shù)的區(qū)間在迭代的過(guò)程中容易發(fā)散，在數(shù)值計(jì)算方面存在困難。目前常見(jiàn)的辦法是將區(qū)間模型修正問(wèn)題轉(zhuǎn)化為確定性修正過(guò)程，如基于區(qū)間響應(yīng)面模型，利用響應(yīng)上下限構(gòu)造目標(biāo)函數(shù)并按照區(qū)間序數(shù)關(guān)系，將不確定性優(yōu)化目標(biāo)和約束條件轉(zhuǎn)化為確定性表達(dá)式，從而避免區(qū)間數(shù)運(yùn)算可能導(dǎo)致的優(yōu)化困難[17]。此外，也有學(xué)者提出將區(qū)間分析與靈敏度分析結(jié)合起來(lái)修正區(qū)間參數(shù)的中值與區(qū)間寬度[18]，但當(dāng)參數(shù)初始區(qū)間估計(jì)不合理時(shí)，二者同時(shí)進(jìn)行修正可能降低優(yōu)化效率，并且可能使修正結(jié)果與實(shí)際結(jié)果出現(xiàn)偏差。

本文考慮試驗(yàn)?zāi)B(tài)數(shù)據(jù)以及結(jié)構(gòu)參數(shù)的區(qū)間不確定性，將不確定性模型修正分成兩步獨(dú)立的優(yōu)化過(guò)程，基于靈敏度分析與響應(yīng)面模型，完成了結(jié)構(gòu)中區(qū)間不確定性參數(shù)進(jìn)行識(shí)別。最后，通過(guò)兩個(gè)算例對(duì)所提出的方法進(jìn)行了驗(yàn)證。

1 試驗(yàn)樣本區(qū)間的中值與半徑分析

考慮結(jié)構(gòu)某些參數(shù)為區(qū)間不確定性時(shí)的動(dòng)力學(xué)模型修正，首先需要確定關(guān)心試驗(yàn)結(jié)果的區(qū)間不確定性，但是，結(jié)構(gòu)動(dòng)力學(xué)試驗(yàn)所獲取的響應(yīng)樣本量通常是小樣本，其樣本數(shù)一般在幾個(gè)～十幾個(gè)。為此，在進(jìn)行模型修正前需要對(duì)小樣本試驗(yàn)結(jié)果的區(qū)間進(jìn)行準(zhǔn)確估計(jì)，以得到可靠性的修正模型。

針對(duì)小樣本數(shù)據(jù)的區(qū)間邊界估計(jì)，常用的方法有Bootstrap法、核密度估計(jì)法、模糊范數(shù)法和灰自助法。有研究表明，變窗寬的核密度估計(jì)方法估計(jì)得到小樣本數(shù)據(jù)的區(qū)間邊界更可靠[19]，因此，本文采用變窗寬的核密度估計(jì)方法來(lái)估計(jì)試驗(yàn)結(jié)果的區(qū)間邊界，以獲得變量區(qū)間中值與半徑。

核密度估計(jì)作為一種非參數(shù)估計(jì)方法，是區(qū)間邊界估計(jì)的一種重要方法。它并不提前假定總體分布的形式，而是以核函數(shù)(如三角形窗、高斯窗、指數(shù)窗等)為權(quán)重，用待估計(jì)點(diǎn)x鄰域內(nèi)其他樣本的值作為該點(diǎn)的估計(jì)值，其估計(jì)函數(shù)如式(1)所示。

(1)

式中：K為核函數(shù)；h為窗寬

(2a)

(2b)

2 不確定性參數(shù)識(shí)別方法

本文基本思想是將不確定性參數(shù)識(shí)別方法分為兩步實(shí)施：第一步利用傳統(tǒng)確定性有限元模型修正與參數(shù)識(shí)別方法，并結(jié)合響應(yīng)面代理模型修正得到待修正參數(shù)的中值；第二步基于代理模型，采用靈敏度分析方法分析得到待修正參數(shù)的半徑，結(jié)合第一步修正結(jié)果，完成結(jié)構(gòu)中區(qū)間參數(shù)的識(shí)別。

2.1 區(qū)間參數(shù)的中值優(yōu)化識(shí)別

該問(wèn)題屬于結(jié)構(gòu)確定性參數(shù)識(shí)別方法，傳統(tǒng)的各種確定性模型修正方法(如神經(jīng)網(wǎng)絡(luò)法、遺傳算法、響應(yīng)面法等)均可使用。考慮到響應(yīng)面法由于具有計(jì)算量小，精度較高，且待修正參數(shù)與響應(yīng)量之間具有顯示函數(shù)關(guān)系，便于后續(xù)進(jìn)一步開(kāi)展基于靈敏度分析的不確定性參數(shù)區(qū)間寬度修正，本文采用多項(xiàng)式響應(yīng)面函數(shù)作為代理模型，開(kāi)展結(jié)構(gòu)系統(tǒng)的不確定性有限元模型修正與參數(shù)識(shí)別。

響應(yīng)面法的基本思想是假設(shè)隨機(jī)輸入變量對(duì)結(jié)構(gòu)響應(yīng)變量的影響可用數(shù)學(xué)函數(shù)來(lái)表達(dá)，通過(guò)確定性有限元方法在隨機(jī)輸入變量空間構(gòu)造有限樣本點(diǎn)，用二次多項(xiàng)式擬合這些樣本點(diǎn)，得到響應(yīng)面函數(shù)。

研究表明，參數(shù)間相互效應(yīng)對(duì)響應(yīng)面模型總方差的貢獻(xiàn)非常小[20]，故本文采用無(wú)交叉項(xiàng)的二次多項(xiàng)式模型，其表達(dá)式為

(3)

式中:b0,bi及bii為待定系數(shù)，需根據(jù)試驗(yàn)樣本點(diǎn)的計(jì)算結(jié)果對(duì)其進(jìn)行識(shí)別。

基于多元線性回歸的多項(xiàng)式響應(yīng)面需數(shù)據(jù)樣本，針對(duì)實(shí)際工程中因素較多的情況，常用的試驗(yàn)設(shè)計(jì)方法有中心復(fù)合設(shè)計(jì)法(又稱CCD(Central Composite Design)法)和Box-Behnken矩陣設(shè)計(jì)法(又稱BBM(Box-Behnken Method)法)。本文采用試驗(yàn)設(shè)計(jì)方法為BBM法，該方法的設(shè)計(jì)點(diǎn)由設(shè)計(jì)參數(shù)的中心點(diǎn)和邊界點(diǎn)組成。

通常響應(yīng)面的精度依據(jù)可根據(jù)式(4)判斷

(4)

傳統(tǒng)確定性的有限元模型修正可歸結(jié)為如下優(yōu)化問(wèn)題

(5)

由于待修正參數(shù)的初始區(qū)間可能估計(jì)不準(zhǔn)，或?qū)τ诟叨确蔷€性模型，由于待修正參數(shù)的初始設(shè)計(jì)空間較大，擬合精度難以保證，參數(shù)識(shí)別可能陷入局部最優(yōu)，因此，借鑒文獻(xiàn)[21]自適應(yīng)響應(yīng)面方法解決該問(wèn)題，即待修正參數(shù)的中值修正可能需要進(jìn)行多步確定性的修正過(guò)程。設(shè)定優(yōu)化收斂的判據(jù)|Fk+1-Fk|/Fk≤ε1或Fk≤ε2，若不滿足，則以最近獲得的參數(shù)優(yōu)化值為中心，縮放或平移新的參數(shù)關(guān)心區(qū)域，重新取得試驗(yàn)設(shè)計(jì)點(diǎn)，并對(duì)結(jié)構(gòu)開(kāi)展相應(yīng)的動(dòng)力學(xué)計(jì)算分析、響應(yīng)面構(gòu)建及優(yōu)化分析。結(jié)構(gòu)待修正參數(shù)的中值修正流程圖，如圖1所示。

圖1 基于響應(yīng)面的區(qū)間參數(shù)中值識(shí)別流程圖Fig.1 Flow chart identifying the mean of uncertain parameters

2.2 區(qū)間參數(shù)的半徑優(yōu)化識(shí)別

假設(shè)待修正參數(shù)p是具有不確定性的結(jié)構(gòu)參數(shù)，用區(qū)間方法可以將其描述為

(6)

對(duì)于具有區(qū)間參數(shù)p的結(jié)構(gòu)，其特征值λ(p)也具有區(qū)間性質(zhì)，由于特征值在前章節(jié)通過(guò)二次響應(yīng)面模型近似，則將其在區(qū)間中值附近進(jìn)行Taylor展開(kāi)，有

λ(p)=λ(pc+Δp)=

(7)

采用區(qū)間數(shù)學(xué)中的區(qū)間自然擴(kuò)張理論，根據(jù)式(7)，通過(guò)區(qū)間運(yùn)算，可得特征值區(qū)間的上、下界及區(qū)間不確定性分別為

(8a)

(8b)

(8c)

采用靈敏度分析的方法迭代求解優(yōu)化問(wèn)題，有限元模型修正問(wèn)題描述為

ye-ya_ j=Sj(pj+1-pj)

(9)

式中：Sj=?ya_ j/?pj為待修正參數(shù)關(guān)于模態(tài)參數(shù)的靈敏度矩陣；可以根據(jù)響應(yīng)面模型式(3)方便計(jì)算得到。

考慮到試驗(yàn)數(shù)據(jù)與結(jié)構(gòu)參數(shù)的區(qū)間不確定性，將式(9)中的變量采用區(qū)間方法表示為

(10)

(11)

(12)

將式(10)代入式(9)，可得

(13)

為了識(shí)別得到待修正參數(shù)的區(qū)間半徑，采用攝動(dòng)法，將式(13)中關(guān)于Δ的零階項(xiàng)和一階項(xiàng)分離，并忽略二階項(xiàng)的影響，得到

(14)

求解式(14)，得到待修正參數(shù)區(qū)間半徑的迭代表達(dá)式為

(15)

根據(jù)上述分析，得到待修正參數(shù)的區(qū)間半徑優(yōu)化識(shí)別過(guò)程，如圖2所示。

圖2 不確定性參數(shù)區(qū)間半徑識(shí)別流程圖Fig.2 Flow chart identifying the radius of uncertain parameters

3 算例研究

3.1 兩自由度系統(tǒng)的參數(shù)識(shí)別

如圖3所示兩自由度系統(tǒng)，模型中確定性的參數(shù)有m1=m2=1 kg，k1=1 N/m。模型中不確定性參數(shù)取值區(qū)間為k2=[0.8,1.2] N/m，k3=[0.9,1.1] N/m。

圖3 兩自由度質(zhì)量—彈簧系統(tǒng)Fig.3 Two degree of freedom mass-spring system

假定不確定性參數(shù)為均勻分布，采用拉丁超立方抽樣方法構(gòu)造20個(gè)試驗(yàn)樣本，計(jì)算得到系統(tǒng)的前兩階固有頻率樣本，并采用非參數(shù)核密度估計(jì)法估計(jì)試驗(yàn)樣本的95%置信區(qū)間，根據(jù)式(2a)、式(2b)計(jì)算得到頻率f1,f2試驗(yàn)結(jié)果的樣本中值分別是0.999 3 rad·s-1和1.731 3 rad·s-1，區(qū)間半徑分別是0.033 2 rad·s-1和0.157 0 rad·s-1。

假定k2,k3的初始區(qū)間值都為[1.5,2.5] N/m。根據(jù)圖1所示區(qū)間參數(shù)中值優(yōu)化流程，開(kāi)展其中值識(shí)別，經(jīng)過(guò)5步迭代即收斂，且每個(gè)迭代步中所建立的響應(yīng)面模型的R2>0.999，k2，k3中值的迭代收斂情況如圖4所示。優(yōu)化識(shí)別得到k2,k3的中值分別為0.999 5 N/m和0.997 3 N/m。

圖4 區(qū)間參數(shù)中值迭代收斂情況Fig.4 Convergence curves of the means of uncertain parameters

根據(jù)圖2所示參數(shù)區(qū)間半徑優(yōu)化流程，通過(guò)3步迭代即收斂，k2,k3半徑的迭代收斂情況如圖5所示。優(yōu)化識(shí)別得到k2,k3的區(qū)間半徑分別是0.239 6 N/m和0.127 6 N/m。

圖6為修正前后不確定性參數(shù)區(qū)間的比較。根據(jù)待修正參數(shù)修正前后的區(qū)間值，采用均勻分布的拉丁超立方抽樣方法構(gòu)造1 000個(gè)隨機(jī)樣本，并計(jì)算系統(tǒng)的前兩階固有頻率，并與試驗(yàn)數(shù)據(jù)進(jìn)行比較，散點(diǎn)圖如圖7所示。

圖5 區(qū)間參數(shù)半徑迭代收斂情況Fig.5 Convergence curves of the radiuses of uncertain parameters

圖6 修正前后區(qū)間參數(shù)的比較Fig.6 Comparison between the original and the updated interval of uncertain parameters

圖7 修正前后系統(tǒng)固有頻率計(jì)算結(jié)果與試驗(yàn)結(jié)果比較Fig.7 Comparison between the simulation and the experiment frequency results

3.2 反射鏡架結(jié)構(gòu)的參數(shù)識(shí)別

神光III光機(jī)裝置是具有48路光束的對(duì)打靶精度有極高要求的慣性約束聚變裝置，該裝置有276個(gè)反射鏡架系統(tǒng)，而反射鏡架在環(huán)境激勵(lì)下的穩(wěn)定性是影響整個(gè)裝置穩(wěn)定性的關(guān)鍵因素[22]。本文考慮如圖8所示的典型鏡架系統(tǒng)，該系統(tǒng)有鏡框組件和鏡架支撐兩部分組成。

鏡架支撐結(jié)構(gòu)的材料為鋼，名義楊氏模量為200 GPa，密度為7 850 kg/m3，泊松比為0.3。鏡框組件中鏡片的材料為K9玻璃，名義楊氏模量為80 GPa，密度為2 510 kg/m3，泊松比為0.21。整個(gè)鏡架系統(tǒng)約1.32 m高。

圖8 反射鏡架結(jié)構(gòu)示意圖Fig.8 Sketch map of the transport mirror system

鏡架系統(tǒng)的不確定性主要是其底部與基礎(chǔ)安裝波動(dòng)以及鏡架支撐焊接工藝的差異帶來(lái)結(jié)構(gòu)剛度的不確定性，導(dǎo)致其動(dòng)力學(xué)特性具有不確定性。為了研究不確定性對(duì)系統(tǒng)固有頻率的影響，對(duì)同類型的鏡架系統(tǒng)開(kāi)展了10次模態(tài)試驗(yàn)。在感興趣的頻率范圍內(nèi)，有三階模態(tài)，其影響結(jié)構(gòu)的振動(dòng)環(huán)境下的動(dòng)力學(xué)響應(yīng)。因此，開(kāi)展動(dòng)力學(xué)響應(yīng)計(jì)算前應(yīng)對(duì)該三階模態(tài)進(jìn)行模型修正。通過(guò)模態(tài)試驗(yàn)獲得結(jié)構(gòu)前兩階彎曲頻率的10個(gè)樣本，如表1所示。

表1 鏡架結(jié)構(gòu)頻率試驗(yàn)值樣本

在ANSYS中采用SOLID185單元及SHELL181單元建立結(jié)構(gòu)的有限元模型，如圖9所示。結(jié)構(gòu)前三階計(jì)算模態(tài)振型如圖10所示。

圖9 鏡架結(jié)構(gòu)的有限元模型Fig.9 FE model of transport mirror system

圖10 鏡架結(jié)構(gòu)前三階模態(tài)振型Fig.10 Modal shapes of transport mirror system

安裝帶來(lái)的剛度不確定性在模型中用鏡架支撐底部材料彈性模量E1表征，焊接不確定性用鏡架支撐底部材料彈性模量E2表征。不確定性參數(shù)識(shí)別就是識(shí)別區(qū)間量E1與E2。

假定E1與E2的初始區(qū)間值分別為[1.0,1.5]×1011Pa 與[1.8,2.0]×1011Pa。首先根據(jù)式(2a)求出試驗(yàn)結(jié)果的中值，再根據(jù)圖1所示區(qū)間參數(shù)中值優(yōu)化流程，開(kāi)展E1與E2的中值識(shí)別，經(jīng)過(guò)6步迭代即收斂，且每個(gè)迭代步中所建立的響應(yīng)面模型的R2>0.999，E1與E2中值迭代收斂情況如圖11所示。

圖11 彈性模量中值收斂情況Fig.11 Convergence curve of the mean of elastic modulus

根據(jù)圖2所示參數(shù)區(qū)間半徑優(yōu)化流程，通過(guò)6步迭代即收斂，E1與E2半徑的迭代收斂情況如圖12所示。

圖12 彈性模量半徑收斂情況圖Fig.12 Convergence curve of the radius of elastic modulus

圖13為修正前后不確定性參數(shù)區(qū)間的比較。根據(jù)待修正參數(shù)修正前后的區(qū)間值，采用均勻分布的拉丁超立方抽樣方法構(gòu)造1 000個(gè)隨機(jī)樣本，并計(jì)算反射鏡架系統(tǒng)的前三階固有頻率，并與試驗(yàn)數(shù)據(jù)進(jìn)行比較，散點(diǎn)圖如圖14所示。并且，根據(jù)1 000個(gè)樣本估算系統(tǒng)前三階固有頻率的95%置信區(qū)間，并將其作為系統(tǒng)前三階固有頻率的區(qū)間，與試驗(yàn)結(jié)果的比較如表2所示。

從圖14與表2可知，修正結(jié)果與試驗(yàn)結(jié)果吻合良好，且修正后的頻率范圍與修正前相比顯得更符合工程實(shí)際。

圖13 不確定性參數(shù)修正前后取值的比較Fig.13 Comparison between the original and the updated interval of uncertain parameters

圖14 修正前后固有頻率計(jì)算結(jié)果與試驗(yàn)結(jié)果比較Fig.14 Scatter map of the first three natural frequencies

模態(tài)試驗(yàn)區(qū)間/Hz初始計(jì)算區(qū)間/Hz誤差/%修正區(qū)間/Hz誤差/%1[21.46,21.80][23.39,24.42][8.99,12.02][21.45,21.79][-0.07,-0.04]2[22.09,22.41][24.05,25.01][8.87,11.60][22.08,22.42][-0.04,0.05]3[41.20,41.81][44.47,46.12][7.94,10.31][40.93,41.55][-0.66,-0.62]

3.2 討 論

有限元模型修正與參數(shù)識(shí)別的時(shí)間主要取決于有限元模型的計(jì)算時(shí)間，在考慮不確定性時(shí)鏡架系統(tǒng)的模型修正僅在第一步中值修正中開(kāi)展了有限元計(jì)算。在構(gòu)建響應(yīng)面的每次迭代步中有9個(gè)樣本，經(jīng)過(guò)6個(gè)迭代步，故總共需要開(kāi)展54次有限元模態(tài)計(jì)算。但是，如果不用代理模型，僅確定性的模型修正就需要開(kāi)展上百次有限元計(jì)算，而考慮不確定性時(shí)，有限元計(jì)算次數(shù)會(huì)更多，因此，本文基于響應(yīng)面與靈敏度分析的區(qū)間不確定性模型修正是高效的。

4 結(jié) 論

(1)不確定性在工程實(shí)際中廣泛存在，針對(duì)考慮不確定性的動(dòng)力學(xué)模型參數(shù)識(shí)別，本文提出了模型中區(qū)間不確定性參數(shù)識(shí)別的分步實(shí)施方法，第一步基于響應(yīng)面模型對(duì)區(qū)間參數(shù)中值進(jìn)行修正，第二步基于響應(yīng)面模型并結(jié)合靈敏度分析實(shí)現(xiàn)區(qū)間參數(shù)的半徑修正，并給出了兩個(gè)修正過(guò)程的實(shí)施流程。由于采用響應(yīng)面模型，可以大幅提高修正效率。

(2)在該思想下開(kāi)展不確定性參數(shù)識(shí)別，在第一步中傳統(tǒng)確定性的動(dòng)力學(xué)參數(shù)識(shí)別與模型修正方法均可使用，便于實(shí)際應(yīng)用。而且，提出首先對(duì)試驗(yàn)小樣本采用核密度估計(jì)方法進(jìn)行合理的區(qū)間估計(jì)，這使得參數(shù)識(shí)別結(jié)果對(duì)工程結(jié)構(gòu)的計(jì)算預(yù)測(cè)更具有普遍指導(dǎo)意義。

(3)針對(duì)本文所提出的區(qū)間不確定性參數(shù)識(shí)別方法，通過(guò)兩個(gè)算例對(duì)參數(shù)識(shí)別方法與流程進(jìn)行了展示，識(shí)別結(jié)果與試驗(yàn)結(jié)果吻合良好。

(4)本文提出的分步實(shí)施的區(qū)間不確定性參數(shù)識(shí)別方法可有效避免區(qū)間優(yōu)化導(dǎo)致的收斂困難，并提高參數(shù)識(shí)別效率，且在實(shí)際工程中易于實(shí)施。