基于分步優(yōu)化的區(qū)間模型修正方法

秦仙蓉 陸慧澄 張曉輝 孫遠(yuǎn)韜 張 氫

同濟(jì)大學(xué)機(jī)械與能源工程學(xué)院 上海 201804

0 引言

工程機(jī)械在國(guó)家現(xiàn)代化發(fā)展和城鎮(zhèn)化建設(shè)中扮演著十分重要的角色。建立結(jié)構(gòu)有限元模型，代替制造周期長(zhǎng)、成本高昂的實(shí)物模型進(jìn)行性能監(jiān)測(cè)以保證其結(jié)構(gòu)安全十分常見[1]。在建立有限元模型的過程中，通常會(huì)對(duì)一些設(shè)計(jì)細(xì)節(jié)進(jìn)行一定的假定和近似處理，這會(huì)造成有限元模型與實(shí)際結(jié)構(gòu)存在差異[2]。因此，有必要進(jìn)行有限元模型修正，使計(jì)算得到的響應(yīng)值與實(shí)測(cè)值趨于一致，以提高精度。

對(duì)于機(jī)械結(jié)構(gòu)而言，由于設(shè)計(jì)制造、裝配及連接方式的不同等原因，會(huì)造成機(jī)械結(jié)構(gòu)的模態(tài)指標(biāo)呈現(xiàn)一定的分散特征，故在工程實(shí)際中考慮不確定性問題對(duì)模型修正有重要意義。若結(jié)構(gòu)參數(shù)為區(qū)間不確定參數(shù)，則機(jī)械結(jié)構(gòu)的模態(tài)頻率也為區(qū)間值，求解這種含區(qū)間參數(shù)的模態(tài)頻率問題即為模態(tài)頻率的區(qū)間分析問題。姜東等[3]在假設(shè)參數(shù)不確定性范圍較小的情況下，推導(dǎo)了參數(shù)中值和參數(shù)半徑的迭代公式，并對(duì)復(fù)合材料板結(jié)構(gòu)進(jìn)行區(qū)間修正。方圣恩等[4]結(jié)合響應(yīng)面模型和區(qū)間算法提出了區(qū)間響應(yīng)面模型，解決了區(qū)間運(yùn)算過程中的區(qū)間擴(kuò)張問題，并通過區(qū)間響應(yīng)的上下界建立優(yōu)化目標(biāo)函數(shù)，避免了復(fù)雜的區(qū)間參數(shù)靈敏度計(jì)算，提高了修正效率。駱勇鵬等[5]利用二次多項(xiàng)式響應(yīng)面和模態(tài)區(qū)間分析理論，提出了區(qū)間逆響應(yīng)面修正方法，修正結(jié)果表明該方法可有效控制區(qū)間過估計(jì)現(xiàn)象，且計(jì)算效率較高。秦仙蓉等[6]基于Kriging替代模型及多目標(biāo)遺傳算法對(duì)岸邊集裝箱起重機(jī)結(jié)構(gòu)進(jìn)行了不確定性有限元模型修正。姜東等[7]基于攝動(dòng)法，針對(duì)參數(shù)不確定性較小的模型修正問題，推導(dǎo)了待修正參數(shù)均值和協(xié)方差矩陣的迭代格式。杜永峰等[8]提出了一種小波Shannon熵與最小標(biāo)準(zhǔn)差相結(jié)合的方法對(duì)復(fù)Morlet小波的中心頻率與帶寬參數(shù)同時(shí)進(jìn)行優(yōu)化，其優(yōu)化后可準(zhǔn)確地識(shí)別結(jié)構(gòu)的密集模態(tài)參數(shù)，并且具有一定的抗噪性能。

本文提出以分步優(yōu)化策略和區(qū)間分析方法為核心的模型修正方法，將區(qū)間修正轉(zhuǎn)化為中值修正與半徑修正分步進(jìn)行，提高修正效率，在半徑修正中采用區(qū)間攝動(dòng)法和Chebyshev多項(xiàng)式法提高修正精度。并將上述方法與替代模型相結(jié)合，解決了該方法應(yīng)用于復(fù)雜有限元模型修正時(shí)重復(fù)調(diào)用有限元模型浪費(fèi)時(shí)間成本的問題。通過三個(gè)算例驗(yàn)證了本文方法的合理性以及精度。

1 不確定性模型修正方法數(shù)學(xué)模型

若有n維區(qū)間參數(shù)向量x={x1，…，xi，…，xn}，根據(jù)區(qū)間分析方法將其分解為中值向量xc與半徑向量Δx和的形式，即

式中：e∈[-1，1]，δx∈[-Δx，Δx]。

本文提出的區(qū)間模型修正方法的核心策略是分步優(yōu)化，先修正參數(shù)中值，再基于修正后的參數(shù)中值，修正參數(shù)半徑，最終得到參數(shù)區(qū)間。該法將區(qū)間修正問題轉(zhuǎn)化為求解參數(shù)中值和參數(shù)半徑的優(yōu)化問題。

參數(shù)中值修正的目標(biāo)函數(shù)為

式中：fc為參數(shù)中值優(yōu)化目標(biāo)函數(shù)，xc為參數(shù)中值，和分別為有限元模型和實(shí)測(cè)的固有頻率中值，m為參與修正的模態(tài)階次，LBc和UBc分別為參數(shù)中值優(yōu)化時(shí)約束優(yōu)化所設(shè)置的上下界。

參數(shù)半徑修正的目標(biāo)函數(shù)為

式中：fr為參數(shù)半徑優(yōu)化目標(biāo)函數(shù)，Δx為參數(shù)半徑，和分別為有限元和實(shí)測(cè)的模態(tài)頻率半徑，m為參與修正的模態(tài)階次，LBr和UBr分別為參數(shù)半徑優(yōu)化時(shí)約束優(yōu)化所設(shè)置的上下界。

結(jié)合圖1所示的基于分步優(yōu)化的區(qū)間模型修正流程圖對(duì)修正過程進(jìn)行描述：1)選取有限元模型修正變量，并構(gòu)建初始區(qū)間；2)通過數(shù)值計(jì)算的頻率響應(yīng)區(qū)間和實(shí)測(cè)的頻率區(qū)間分別構(gòu)建參數(shù)區(qū)間中值和半徑的優(yōu)化目標(biāo)函數(shù)；3)通過一次確定性修正對(duì)參數(shù)中值進(jìn)行修正優(yōu)化；4)通過2種區(qū)間分析方法對(duì)參數(shù)半徑進(jìn)行不確定性修正優(yōu)化；5)得到參數(shù)區(qū)間，完成修正。

圖1 基于分步優(yōu)化的區(qū)間模型修正流程圖

2 區(qū)間分析方法

在修正參數(shù)半徑時(shí)，本文引入基于區(qū)間攝動(dòng)法和Chebyshev多項(xiàng)式的區(qū)間分析方法對(duì)式(3)所示的目標(biāo)函數(shù)進(jìn)行求解，其中區(qū)間攝動(dòng)法能有效地提高修正效率，Chebyshev多項(xiàng)式在區(qū)間范圍內(nèi)具有最佳平方逼近的性質(zhì)，通過在區(qū)間范圍內(nèi)進(jìn)行插值計(jì)算可近似得到輸出的區(qū)間范圍，且精度高[9]。

2.1 基于區(qū)間攝動(dòng)法的區(qū)間分析法

區(qū)間攝動(dòng)法[10]的核心是函數(shù)一階泰勒展開，它僅需計(jì)算響應(yīng)在某一點(diǎn)處偏導(dǎo)數(shù)矩陣，無(wú)需進(jìn)行多次有限元計(jì)算來(lái)判斷頻率的上下界。

將響應(yīng)f在區(qū)間參數(shù)中值xc處進(jìn)行泰勒級(jí)數(shù)展開，若x的變化區(qū)間較小，僅保留其常數(shù)項(xiàng)和一階項(xiàng)[11]，如式(4)所示

對(duì)于式(4)中的偏導(dǎo)數(shù)矩陣求解問題，常需要先行求解剛度矩陣和質(zhì)量矩陣的偏導(dǎo)數(shù)矩陣。然而對(duì)于復(fù)雜的機(jī)械結(jié)構(gòu)而言，其剛度和質(zhì)量矩陣無(wú)法得到[12,13]。本文利用中心差分法來(lái)近似求解響應(yīng)對(duì)參數(shù)的偏導(dǎo)數(shù)矩陣。可以得到響應(yīng)區(qū)間的半徑為

2.2 基于Chebyshev多項(xiàng)式的區(qū)間分析法

Chebyshev多項(xiàng)式[14]無(wú)需知道參數(shù)與頻率的單調(diào)性關(guān)系，可通過求解張量積得到Chebyshev多項(xiàng)式系數(shù)，再基于插值點(diǎn)和多項(xiàng)式公式得到頻率的區(qū)間范圍。

對(duì)于區(qū)間變量xI，則f(x)的Chebyshev多項(xiàng)式區(qū)間表達(dá)為

式中：k為Chebyshev多項(xiàng)式階次，θI=[0，π]。

因此，函數(shù)f(x)的求解問題轉(zhuǎn)換成求解其Chebyshev多項(xiàng)式系數(shù)fi的問題，根據(jù)Chebyshev基函數(shù)的正交性和Gaussie-Chebyshev積分公式可得

式中：xj為Chebyshev多項(xiàng)式的插值點(diǎn)，xj=cosθj(j=1，2，…，p)；p為插值點(diǎn)數(shù)量，p=k+1。

上述將函數(shù)f(x)展開成Chebyshev多項(xiàng)式是一維問題，即只有一個(gè)區(qū)間不確定性參數(shù)。對(duì)于有n個(gè)區(qū)間不確定性參數(shù)而言，其Chebyshev多項(xiàng)式可利用n個(gè)一維Chebyshev多項(xiàng)式的張量積形式來(lái)構(gòu)造，即

式中：θi=arccosxi∈ [0，π]。

則n維k階Chebyshev多項(xiàng)式的公式為

根據(jù)式(6)將式(9)轉(zhuǎn)化成區(qū)間函數(shù)形式

式中：fi1,…,in為Chebyshev多項(xiàng)式系數(shù)，l為Ci1,i2,…,ik下標(biāo)中含0的個(gè)數(shù)。

可得響應(yīng)區(qū)間的半徑

3 均勻簡(jiǎn)支梁數(shù)值算例

本節(jié)采用一根均勻簡(jiǎn)支梁[15]作為數(shù)值算例，如圖2所示。并與其他文獻(xiàn)中的區(qū)間模型修正方法進(jìn)行了對(duì)比。

圖2 均勻簡(jiǎn)支梁結(jié)構(gòu)

均勻簡(jiǎn)支梁的各參數(shù)的取值如表1所示。

表1 均勻簡(jiǎn)支梁結(jié)構(gòu)參數(shù)

根據(jù)文獻(xiàn)[15]的靈敏度分析結(jié)果，選取梁的密度和截面寬度作為不確定性待修正參數(shù)。假設(shè)兩個(gè)參數(shù)的真實(shí)區(qū)間為[2 125, 2 375]和[0.263,0.288]，利用拉丁超立方抽樣方法在真實(shí)區(qū)間內(nèi)抽樣50次并計(jì)算得到真實(shí)頻率區(qū)間，如表2所示。

表2 均勻簡(jiǎn)支梁模態(tài)頻率的真實(shí)區(qū)間 Hz

兩個(gè)參數(shù)的初始區(qū)間設(shè)為[2 000，3 000]、[0.2，0.3]，根據(jù)圖1步驟對(duì)參數(shù)中值和參數(shù)半徑進(jìn)行修正，得到修正后的參數(shù)區(qū)間。將修正結(jié)果與區(qū)間響應(yīng)面模型直接修正參數(shù)區(qū)間上下限的[4]修正結(jié)果進(jìn)行對(duì)比，如表3所示。

表3 均勻簡(jiǎn)支梁修正后的的參數(shù)區(qū)間

區(qū)間響應(yīng)面擬合流程如圖3所示，參數(shù)為密度和截面寬度，參數(shù)區(qū)間為初始區(qū)間，通過拉丁超立方抽樣100次構(gòu)建響應(yīng)面。

圖3 區(qū)間響應(yīng)面構(gòu)建方法

通過對(duì)比，可以看出：本文方法比區(qū)間響應(yīng)面模型修正方法精度更高。經(jīng)過本文提出的方法修正后，區(qū)間上下界最大誤差由26.3%降低到了1.0%以下。

4 工程案例

本文以圖4所示的某型塔式起重機(jī)（以下簡(jiǎn)稱塔機(jī)）為對(duì)象構(gòu)建有限元模型。塔身、起重臂、平衡臂、塔頂?shù)冉饘俳Y(jié)構(gòu)采用梁?jiǎn)卧M(jìn)行建模；起重臂和平衡臂的拉桿采用桿單元進(jìn)行建模；將標(biāo)準(zhǔn)節(jié)之間與加強(qiáng)節(jié)之間采用固結(jié)處理，起重臂各節(jié)臂之間，拉桿與塔頂之間、兩臂與回轉(zhuǎn)塔身之間采用鉸接連接。將附屬構(gòu)件以等效質(zhì)量的形式設(shè)置在相應(yīng)的位置上。

圖4 某型塔機(jī)

在修正前需選定不確定參數(shù)，基于塔機(jī)有限元模型，選取塔身標(biāo)準(zhǔn)節(jié)主肢長(zhǎng)度、平衡臂平臺(tái)端相互鉸接的斜腹底邊距離、上塔身水平腹桿長(zhǎng)度、塔身水平腹桿長(zhǎng)度、起重臂下弦桿壁厚、平衡臂尾部主梁高度、塔身主肢截面寬度、塔身主肢截面厚度、彈性模量、泊松比，10個(gè)參數(shù)進(jìn)行Sobol全局靈敏度分析，對(duì)應(yīng)前5階頻率每個(gè)參數(shù)都有5個(gè)靈敏度指標(biāo)，其綜合靈敏度采用加權(quán)疊加的形式計(jì)算，權(quán)重系數(shù)皆取0.2，靈敏度分析結(jié)果如圖5所示。選取4個(gè)靈敏度較高的參數(shù)作為不確定性待修正參數(shù)，以有限元模型各參數(shù)初始值為中值，以4個(gè)參數(shù)的初始中值為基準(zhǔn)，對(duì)4個(gè)參數(shù)歸一化進(jìn)行歸一化處理后，皆取[0.75, 1.25]為4個(gè)參數(shù)的初始區(qū)間，如表4所示。

表4 塔機(jī)結(jié)構(gòu)不確定性修正參數(shù)

圖5 塔機(jī)修正變量對(duì)前5階頻率的靈敏度

與第3小節(jié)梁算例不同的是：塔機(jī)是復(fù)雜機(jī)械結(jié)構(gòu)，其響應(yīng)頻率區(qū)間求解依賴于有限元模型，如在修正時(shí)反復(fù)調(diào)用有限元模型進(jìn)行目標(biāo)函數(shù)計(jì)算顯然是效率低下的，尤其是Chebyshev多項(xiàng)式法需要計(jì)算參數(shù)插值點(diǎn)組合時(shí)的頻率，計(jì)算成本呈指數(shù)級(jí)上升。

針對(duì)上述問題，在塔機(jī)案例中將本文提出的區(qū)間修正方法和替代模型相結(jié)合進(jìn)行修正，替代模型選擇梁算例中用于對(duì)比的區(qū)間響應(yīng)面模型。引入?yún)^(qū)間響應(yīng)面模型后的修正步驟與圖1無(wú)異，僅在求解中值與半徑目標(biāo)函數(shù)時(shí)，以區(qū)間響應(yīng)面模型替代有限元模型完成求解。同樣將修正結(jié)果與區(qū)間響應(yīng)面直接修正區(qū)間上下限的方法[4]進(jìn)行對(duì)比。

區(qū)間響應(yīng)面的構(gòu)建方法如下，通過拉丁超立方抽樣方法在4個(gè)參數(shù)初始區(qū)間內(nèi)抽樣100次，調(diào)用有限元模型計(jì)算樣本對(duì)應(yīng)的頻率，擬合各階頻率的區(qū)間響應(yīng)面模型，該模型后續(xù)將同時(shí)運(yùn)用于3種修正方法中。

4.1 塔機(jī)實(shí)例

在該型塔機(jī)工作狀態(tài)下，采用無(wú)線環(huán)境激勵(lì)試驗(yàn)?zāi)B(tài)測(cè)試分析系統(tǒng)DH5907A對(duì)其進(jìn)行風(fēng)振響應(yīng)實(shí)測(cè)，得到塔機(jī)實(shí)測(cè)前5階模態(tài)頻率中值與半徑如表5所示。

表5 塔機(jī)實(shí)測(cè)頻率中值和半徑（實(shí)例） Hz

結(jié)合表5中的實(shí)測(cè)頻率中值和頻率中值的響應(yīng)面模型建立中值與半徑修正的目標(biāo)函數(shù)，按圖1所示的流程完成參數(shù)中值與半徑的修正。再用區(qū)間響應(yīng)面模型直接修正區(qū)間上下限[4]，3種方法修正后的參數(shù)區(qū)間如表6所示。

由表6可知，參數(shù)2經(jīng)過本文方法修正后的區(qū)間為點(diǎn)區(qū)間（即區(qū)間半徑為零，僅有區(qū)間中值），其原因是實(shí)測(cè)頻率的波動(dòng)非常小。

表6 塔機(jī)結(jié)構(gòu)修正后的的參數(shù)區(qū)間（實(shí)例）

由于塔機(jī)實(shí)例中的已知條件為真實(shí)頻率區(qū)間而非真實(shí)參數(shù)區(qū)間，為了驗(yàn)證修正效果，通過參數(shù)區(qū)間計(jì)算對(duì)應(yīng)的頻率區(qū)間，并轉(zhuǎn)化成中值與區(qū)間半徑的表達(dá)形式如表7所示。

表7 塔機(jī)結(jié)構(gòu)修正后的頻率中值和頻率半徑（實(shí)例） Hz

由表7可知，3種方法修正的頻率區(qū)間中值相同，這是實(shí)測(cè)頻率的波動(dòng)非常小，模型修正的不確定性也非常小，導(dǎo)致最終修正后得到的頻率中值非常接近；頻率區(qū)間的半徑由于實(shí)測(cè)頻率波動(dòng)過小，3種方法修正效果皆不理想，但可以明顯看出本文提出的方法結(jié)合區(qū)間響應(yīng)面模型后的修正結(jié)果體現(xiàn)了第四階實(shí)測(cè)頻率波動(dòng)最大的趨勢(shì)，而區(qū)間模型直接修正上下限則無(wú)法做到，且除第三階實(shí)測(cè)頻率外，修正精度都更高。

4.2 塔機(jī)數(shù)值算例

由于塔機(jī)的實(shí)測(cè)頻率波動(dòng)程度較小，前述雖體現(xiàn)了本文提出的以分步優(yōu)化策略和區(qū)間分析方法為核心的不確定性模型修正方法結(jié)合替代模型的可行性，但還無(wú)法體現(xiàn)其在修正復(fù)雜結(jié)構(gòu)區(qū)間模型時(shí)的精度。

因此，根據(jù)梁算例的流程構(gòu)造一個(gè)塔機(jī)數(shù)值算例以進(jìn)一步驗(yàn)證本文方法的精度與適用性。為保證算例具備可信度，不確定性參數(shù)仍如表4所示，以4個(gè)參數(shù)的初始中值為基準(zhǔn)，對(duì)3個(gè)參數(shù)歸一化進(jìn)行歸一化處理后，分別取 [0.94, 1]、[0.94, 1]、[0.94, 1]、[0.97, 1.03]為 4個(gè)參數(shù)的真實(shí)區(qū)間，利用拉丁超立方抽樣方法在參數(shù)的真實(shí)區(qū)間內(nèi)抽樣50次得到塔機(jī)前5階頻率的中值和半徑表8所示。

表8 塔機(jī)實(shí)測(cè)頻率中值和半徑（數(shù)值算例）Hz

按照?qǐng)D1中的修正流程，以與4.1小節(jié)同樣的3種方法數(shù)值算例中的4個(gè)不確定性參數(shù)進(jìn)行修正，所得的修正后塔機(jī)結(jié)構(gòu)的參數(shù)區(qū)間如表9所示，本文提出的方法結(jié)合區(qū)間響應(yīng)面模型后的修正精度比以相同的區(qū)間響應(yīng)面模型直接修正區(qū)間上下限更高，其誤差都在0.5%以下，而反觀區(qū)間響應(yīng)面模型直接修正區(qū)間上下限的方法，第1個(gè)參數(shù)與第3個(gè)參數(shù)的修正誤差都大于1.0%。

表9 塔機(jī)結(jié)構(gòu)修正后的的參數(shù)區(qū)間（數(shù)值算例）

5 結(jié)論

本文提出一種以分步優(yōu)化策略和區(qū)間分析方法為核心的模型修正方法，根據(jù)分步優(yōu)化策略，首先修正參數(shù)中值，其次用區(qū)間攝動(dòng)法與Chebyshev多項(xiàng)式法修正參數(shù)半徑，完成區(qū)間模型修正。得到如下結(jié)論：

1）本文所提出的區(qū)間模型修正方法可高精度地修正含區(qū)間參數(shù)的均勻簡(jiǎn)支梁模型，且修正效果優(yōu)于區(qū)間響應(yīng)面模型直接修正參數(shù)區(qū)間上下限的修正方法。

2）本文提出的區(qū)間模型修正方法可應(yīng)用于復(fù)雜結(jié)構(gòu)的模型修正中，為避免反復(fù)調(diào)用有限元模型浪費(fèi)計(jì)算成本，可將本文方法與替代模型結(jié)合使用，運(yùn)用替代模型代為計(jì)算每一步迭代中的目標(biāo)函數(shù)。由本文兩個(gè)塔機(jī)算例可看出，結(jié)合本文方法與區(qū)間響應(yīng)面模型極大提升了修正塔機(jī)模型的效率，且精度高于以同樣的區(qū)間響應(yīng)面模型直接修正區(qū)間上下限的方法。在實(shí)測(cè)頻率響應(yīng)區(qū)間較小的情況下，其修正結(jié)果仍能反映頻率的真實(shí)波動(dòng)趨勢(shì)。