任銘澤，鄧忠民，國兆普

(1.北京航空航天大學(xué) 宇航學(xué)院，北京 100191；2.中國航天科工集團(tuán) 北京動(dòng)力機(jī)械研究所，北京 100074)

近年來，有限元模型修正方法被廣泛應(yīng)用于航空、航天、土木建筑等領(lǐng)域[1]，這類方法利用測(cè)量數(shù)據(jù)作為參考數(shù)據(jù)來提高模型精度，通過識(shí)別和改進(jìn)建模參數(shù)，從而使模型能更好地對(duì)結(jié)構(gòu)的響應(yīng)進(jìn)行預(yù)測(cè)[2]。

確定性結(jié)構(gòu)有限元模型修正方法[3-4]目前已被廣泛應(yīng)用于工程結(jié)構(gòu)中，但在實(shí)際工程結(jié)構(gòu)中普遍存在不確定性[5]。工程中的不確定性主要分為三類[6]：①結(jié)構(gòu)參數(shù)不確定性，例如某些復(fù)合材料的物理參數(shù)存在明顯不確定，螺栓、鉸接等連接結(jié)構(gòu)的邊界條件難以精準(zhǔn)化等；②分析模型不確定性，主要是由仿真軟件本身的邊界條件難以精準(zhǔn)化等；③試驗(yàn)測(cè)量不確定性,例如在測(cè)量過程中，由于人為或者儀器本身誤差導(dǎo)致的不確定性。結(jié)構(gòu)不確定性在很大程度上影響動(dòng)力學(xué)分析結(jié)果的可信度，因此開展考慮不確定性的結(jié)構(gòu)有限元模型修正的研究具有極大的工程意義。

不確定性研究方法主要分為概率與非概率方法[7]。概率方法通常需要大量的試驗(yàn)數(shù)據(jù)來獲得較準(zhǔn)確的統(tǒng)計(jì)信息，但是當(dāng)試驗(yàn)數(shù)據(jù)樣本數(shù)量較少，不足以準(zhǔn)確的描述結(jié)構(gòu)的統(tǒng)計(jì)特性時(shí)，會(huì)引起較大的誤差[8]。非概率方法中具有代表性的就是區(qū)間方法，僅需要參數(shù)的上下界即可進(jìn)行相關(guān)研究，對(duì)于處理少量試驗(yàn)數(shù)據(jù)或某些特殊的不確定性問題有很好的效果。得益于區(qū)間理論的發(fā)展以及Chen等[9]、Ben-Haim等[10]學(xué)者的研究，區(qū)間方法已成功應(yīng)用于結(jié)構(gòu)靜力學(xué)、動(dòng)力學(xué)分析問題。

在不確定性有限元模型修正研究領(lǐng)域，Collins等[11]首先提出的最小方差法；之后，Beck等[12]提出了貝葉斯統(tǒng)計(jì)框架，利用動(dòng)態(tài)響應(yīng)數(shù)據(jù)定量分析結(jié)構(gòu)模型的預(yù)測(cè)精；王登剛等[13]提出了一種結(jié)構(gòu)計(jì)算區(qū)間反演方法，將模型修正問題轉(zhuǎn)化為全局優(yōu)化問題；姜東等[14]提出了在區(qū)間參數(shù)結(jié)構(gòu)特征值分析理論基礎(chǔ)上，利用參數(shù)靈敏度分析進(jìn)行區(qū)間模型修正方法。

線性結(jié)構(gòu)的不確定問題研究已經(jīng)取得很大成果，非線性結(jié)構(gòu)與線性結(jié)構(gòu)相比，非線性結(jié)構(gòu)的響應(yīng)與輸入之間的傳遞函數(shù)會(huì)隨著激勵(lì)的變化而變化，因此非線性響應(yīng)的求解往往會(huì)有更大的計(jì)算量甚至?xí)惺諗啃詥栴}。在不確定非線性分析研究上，張義民等[15]提出了一種概率攝動(dòng)有限元法，并用于解決多自由度非線性隨機(jī)參數(shù)振動(dòng)響應(yīng)分析，但是此方法對(duì)參數(shù)先驗(yàn)信息要求較高；邱志平等[16]采用一種以區(qū)間數(shù)學(xué)為基礎(chǔ),將不確定性參數(shù)用區(qū)間進(jìn)行定量化,借助一階泰勒級(jí)數(shù),給出了近似估計(jì)非線性振動(dòng)系統(tǒng)動(dòng)力響應(yīng)范圍的區(qū)間分析方法，相較于概率攝動(dòng)法結(jié)果略顯保守。已有的非線性動(dòng)力學(xué)研究主要集中在研究已知不確定參數(shù)的區(qū)間分布信息，對(duì)參數(shù)的響應(yīng)區(qū)間的求解問題上，對(duì)不確定非線性模型修正問題還缺少相關(guān)的研究。

本文介紹了一種區(qū)間模型修正法，它是一階區(qū)間攝動(dòng)技術(shù)與傳統(tǒng)的有限元模型修正中代價(jià)函數(shù)相結(jié)合的產(chǎn)物。傳統(tǒng)的有限元模型修正問題將仿真響應(yīng)與實(shí)際響應(yīng)的差的平方和作為優(yōu)化的目標(biāo)函數(shù)，本文基于傳統(tǒng)模型修正方法，結(jié)合區(qū)間攝動(dòng)理論，進(jìn)行不確定性模型修正問題中非概率區(qū)間模型修正問題研究，相較于以往不確定非線性分析研究，降低了對(duì)參數(shù)區(qū)間分布的先驗(yàn)信息需求的同時(shí)，研究通過區(qū)間攝動(dòng)法求解響應(yīng)區(qū)間，并在此基礎(chǔ)上對(duì)不確定參數(shù)區(qū)間進(jìn)行修正。

1 相關(guān)理論

1.1 非線性結(jié)構(gòu)動(dòng)力學(xué)區(qū)間攝動(dòng)理論

考慮N自由度非線性結(jié)構(gòu)系統(tǒng)動(dòng)力學(xué)微分方程

(1)

式中，B為不確定性參數(shù)，可以包含不確定載荷參數(shù)和不確定結(jié)構(gòu)參數(shù)。

(3)

式中，Kc為非線性剛度系數(shù)，用于描述非線性恢復(fù)力，Kc通常用非線性系數(shù)c乘以結(jié)構(gòu)不確定剛度K來表示。

采用區(qū)間攝動(dòng)理論將其分解為區(qū)間均值和區(qū)間半徑兩部分

K=[Kc-ΔK,Kc+ΔK]

(4)

式中：Kc為區(qū)間均值；ΔK為區(qū)間的半徑。

對(duì)含不確定參數(shù)的非線性結(jié)構(gòu)時(shí)域響應(yīng)，相應(yīng)地有

X(t)∈X(t)I=[X0-ΔX,X0+ΔX]

(5)

對(duì)于自由度為N的系統(tǒng)有X=[X1,X2,…，Xn]

X的區(qū)間均值為

(6)

式中，s為樣本點(diǎn)數(shù)量。

區(qū)間半徑為

(7)

時(shí)域動(dòng)響應(yīng)分解為

(8)

1.2 基于區(qū)間攝動(dòng)的模型修正理論

定義目標(biāo)函數(shù)

minr(B)=(Xi(B)-Xi)2

(9)

式中：Xi(B)為仿真值;Xi為真實(shí)值。將不確定參數(shù)定義為確定部分和不確定小量的和的形式

B=Bc+δB

(10)

與此同時(shí)將響應(yīng)也分為確定參數(shù)和不確定參數(shù)兩部分的響應(yīng)

(11)

相應(yīng)的將優(yōu)化的目標(biāo)函數(shù)分為由確定參數(shù)(區(qū)間均值)和不確定參數(shù)(區(qū)間半徑)導(dǎo)致的響應(yīng)兩部分的和的形式

r(B)=r0(Bc)+δr(δB)

(12)

將式(8)、式(11)、式(12)代入式(9)

(13)

式(13)展開并分離為確定參數(shù)部分零階式的式(14)和不確定參數(shù)部分一階式的式(15)

(14)

(15)

(16)

通過求解式(14)中的極小化問題，可以得到不確定參數(shù)均值的估計(jì)?，F(xiàn)在的主要目標(biāo)是找到求解參數(shù)的“真”區(qū)間半徑，使結(jié)構(gòu)響應(yīng)的修正值與試驗(yàn)值相匹配，用ΔXi近似替換式(16)中的δXi作為真實(shí)響應(yīng)的區(qū)間半徑，得到

(17)

通過建立和求解式(14)和式(17)的最小值問題?？梢缘玫讲淮_定參數(shù)的修正區(qū)間，并用均值和區(qū)間半徑區(qū)間來表示。通過上述分析，區(qū)間有限元模型修正過程主要包括兩個(gè)步驟：①區(qū)間均值優(yōu)化；②在均值修正基礎(chǔ)上，區(qū)間半徑優(yōu)化。

2 模型修正算法

2.1 區(qū)間均值修正

區(qū)間均值的修正問題，與確定性模型修正問題類似，可分為以下幾步：

步驟2代入結(jié)構(gòu)動(dòng)力學(xué)時(shí)域動(dòng)響應(yīng)計(jì)算的算法中，得到時(shí)域響應(yīng)的區(qū)間均值的X0(θc)，建立參數(shù)區(qū)間均值和響應(yīng)區(qū)間均值的函數(shù)關(guān)系

(18)

步驟3利用式(14)得到待修正參數(shù)的確定部分(區(qū)間均值)導(dǎo)致的響應(yīng)的與試驗(yàn)數(shù)據(jù)的殘差

ε1=X(θ0)-Xm

(19)

步驟4將均值的修正轉(zhuǎn)化為了如下的優(yōu)化問題

(20)

式中,下標(biāo)l和u分別為修正參數(shù)的空間的下限和上限，通過對(duì)式(20)進(jìn)行最小化尋優(yōu)可得到參數(shù)的均值的修正值

θc=[θ1c,θ2c,…,θnc]T

(21)

2.2 區(qū)間半徑修正

區(qū)間半徑的修正算法主要有以下幾個(gè)步驟：

步驟1在區(qū)間均值修正的基礎(chǔ)上，基于對(duì)待修正參數(shù)區(qū)間采用固定抽樣間隔的均勻抽樣方法得到樣本點(diǎn)，代入有限元模型計(jì)算得到時(shí)域響應(yīng)區(qū)間；

步驟2建立參數(shù)區(qū)間半徑和響應(yīng)區(qū)間半徑的函數(shù)關(guān)系

δX=δX(δθ)=[δx1(δθ),δx2(δθ),…,δxn(δθ)]T

(22)

其中,

(23)

(24)

步驟3利用式(17)得到待修正參數(shù)不確定部分(區(qū)間半徑)導(dǎo)致的響應(yīng)的與試驗(yàn)數(shù)據(jù)區(qū)間半徑的殘差

ε2=δX(δθ)-ΔXm

(25)

步驟4因此將區(qū)間半徑的修正轉(zhuǎn)化為了如下的優(yōu)化問題

(26)

步驟5通過對(duì)式(26)進(jìn)行最小化尋優(yōu)可得到參數(shù)的區(qū)間半徑的修正值

Δθ=[δθ1,δθ2,…,δθn]T

(27)

步驟6因此待修正參數(shù)區(qū)間的修正值為

Θ=[θc-Δθ,θc+Δθ]

(28)

具體流程如圖1所示。其中區(qū)間均值修正步驟中，均值附近的固定取樣區(qū)間設(shè)定為真實(shí)區(qū)間半徑，灰色部分為真實(shí)的區(qū)間信息。可以看出算法分為兩步：①求解區(qū)間均值修正值；②在此基礎(chǔ)上，修正攝動(dòng)量(區(qū)間半徑)，這保證區(qū)間半徑是在區(qū)間均值修正基礎(chǔ)上的的修正結(jié)果，這是兩個(gè)相關(guān)聯(lián)的變量。

圖1 基于攝動(dòng)法的區(qū)間模型修正流程

3 算例分析

3.1 算例一 兩自由度剛架結(jié)構(gòu)

算例模型采用如圖2所示的兩層剛架結(jié)構(gòu)，鋼架結(jié)構(gòu)的支柱只提供彎曲剛度不計(jì)質(zhì)量，水平桿具有質(zhì)量，只發(fā)生水平位移，系統(tǒng)的運(yùn)動(dòng)微分方程為

圖2 兩自由度剛架結(jié)構(gòu)

(29)

其中,

Kc=c×K×104

(30)

確定的結(jié)構(gòu)參數(shù)包括剛架m1,m2質(zhì)量分別為1 kg，2 kg，等效非線性系數(shù)為Kc，其中c分別取0和1，對(duì)比分析線性和非線性模型修正方法的修正效果，外作用力F(t)=10·sin(100·t)。

假設(shè)不確定的參數(shù)真實(shí)值區(qū)間為

K1=K2=[0.8,1.2]·103N/m。

區(qū)間均勻抽樣的取樣間隔r分別設(shè)定為0.01和0.001，樣本點(diǎn)個(gè)數(shù)由區(qū)間長度除以取樣間隔r得到。時(shí)域動(dòng)響應(yīng)的計(jì)算采用Newmark-β直接積分法，仿真時(shí)間為1 s，仿真步長為0.001 s。

分別假設(shè)待修正的參數(shù)的初始區(qū)間為

initial-1:K1=K2=0.9+[-0.4,0.4]·103N/m

initial-2:K1=K2=0.8+[-0.5,0.5]·103N/m

當(dāng)非線性系數(shù)取0和1時(shí)，區(qū)間的均值和區(qū)間半徑的迭代收斂曲線，如圖3和圖4所示。從圖可知，區(qū)間均值和半徑經(jīng)過二十余次迭代后收斂到真實(shí)值附近，而且當(dāng)初始區(qū)間不同的情況下，修正后參數(shù)區(qū)間均值和半徑均能收斂到真實(shí)區(qū)間附近，這表明在非線性和線性不確定模型修正問題上，該算法具有很好的收斂效果，且算法穩(wěn)定性好。c取0和1兩種情況下修正前后參數(shù)的區(qū)間誤差，如表1、表2和表3所示。從表中數(shù)據(jù)可知，兩種結(jié)構(gòu)形式下，修正后區(qū)間均值和區(qū)間半徑的誤差小于3%，較初始誤差有大幅度下降，說明該方法對(duì)線性和非線性結(jié)構(gòu)均有很高修正精度，表2和表3分別為不同取樣間隔r為0.01和0.001的條件下，區(qū)間修正誤差對(duì)比，結(jié)果表明取樣間隔越小，修正后參數(shù)區(qū)間誤差更小。

表1 修正前后不確定參數(shù)區(qū)間誤差(c=0)

表2 修正前后不確定參數(shù)區(qū)間誤差(c=1，initial-1,r=0.01)

表3 修正前后不確定參數(shù)區(qū)間誤差(c=1,initial-1,r=0.001)

圖3 區(qū)間均值收斂過程

圖4 區(qū)間半徑收斂過程

對(duì)修正前后的參數(shù)區(qū)間，通過均勻取樣的方法構(gòu)造參數(shù)樣本，帶入確定的有限元模型計(jì)算得到修正前后的時(shí)域響應(yīng)的區(qū)間分布。圖5為(c=1，initial-1,r=0.001)條件下，修正前后自由度X2的時(shí)域響應(yīng)的響應(yīng)區(qū)間與真實(shí)的響應(yīng)區(qū)間的對(duì)比，其中:U和 L分別為時(shí)域響應(yīng)區(qū)間的上界和下界；R，Q，H分別為真實(shí)值、修正前和修正后。從圖5可知，修正后的響應(yīng)區(qū)間和真實(shí)的響應(yīng)區(qū)間的曲線基本重合，兩個(gè)區(qū)間在不同采樣時(shí)間都能很好地吻合，能夠準(zhǔn)確反映仿真試驗(yàn)結(jié)果。

圖5 修正前后時(shí)域響應(yīng)X2(C=1, initial-1,r=0.001)

3.2 算例二 非線性懸臂梁

算例模型采用如圖6所示的自由端帶有非線性彈簧的懸臂梁結(jié)構(gòu)，梁的截面為圓形,梁的結(jié)構(gòu)參數(shù)如表4所示，以梁的材料的彈性模量E和材料密度ρ作為不確定參數(shù)的待修正參數(shù)。

圖6 懸臂梁有限元模型

表4 梁結(jié)構(gòu)參數(shù)

假設(shè)梁的材料特性參數(shù)中彈性模量E和材料密度ρ的真實(shí)值分別為

E=0.7+[-0.1,0.1](102GPa)

ρ=2.7+[-0.3,0.3](103kg/m3)

假設(shè)彈性模量E和材料密度ρ的初始值分別為

E=0.49+[-0.2,0.2](102GPa)

ρ=2.43+[-0.4,0.4](103kg/m3)

即彈性模量E的區(qū)間均值被低估了30%，區(qū)間半徑初始誤差+100%；材料密度ρ的區(qū)間均值被低估了10%，區(qū)間半徑初始誤差33.3%。

取立方非線性和分段非線性兩種非線性形式進(jìn)行仿真，對(duì)于立方非線性，非線性彈簧剛度系數(shù)Kc定義為k·c，c為非線性系數(shù)定義為1×103，彈簧恢復(fù)力位移呈現(xiàn)三次關(guān)系，結(jié)構(gòu)的整體運(yùn)動(dòng)學(xué)方程為

(31)

對(duì)于分段非線性，取分段間隔|δ|=1×10-4，非線性恢復(fù)力呈現(xiàn)分段性，非線性恢復(fù)力為

(32)

結(jié)構(gòu)運(yùn)動(dòng)學(xué)方程為

(33)

外力F(t)=600sin(400×t)，梁單元為歐拉-伯努利梁單元。采用Newmark-β法進(jìn)行時(shí)域動(dòng)響應(yīng)的計(jì)算，取仿真時(shí)間t=5 s，時(shí)間步長為0.001 s。采用間隔r為0.001的均勻取樣方法分別對(duì)初始和真實(shí)參數(shù)區(qū)間進(jìn)行抽樣，構(gòu)造參數(shù)樣本空間。

修正后不確定參數(shù)的區(qū)間誤差比較，如表5和表6所示。非線性形式為立方非線性和分段非線性時(shí)，修正后參數(shù)的區(qū)間均值誤差均下降到小于0.1%，彈性模量E區(qū)間半徑誤差從100%分別下降到1.4%和2.0%，材料密度ρ的區(qū)間半徑誤差下降到0.6%和0.8%。兩種非線性形式下不確定參數(shù)的修正后誤差與修正前相比均有大幅的下降。

表5 修正后參數(shù)的區(qū)間誤差(立方非線性,r=0.001)

表6 修正后參數(shù)的區(qū)間誤差(分段非線性,r=0.001)

彈性模量E和材料密度ρ的區(qū)間均值和區(qū)間半徑的迭代收斂曲線(取樣間隔r為0.001)，如圖7和圖8所示。從圖中可以看出兩種非線性形式的兩個(gè)不確定參數(shù)的區(qū)間均值和區(qū)間半徑在二十多次迭代后收斂于真實(shí)值附近。

圖7 彈性模量E收斂過程：區(qū)間均值，區(qū)間半徑

圖8 材料密度ρ收斂過程：區(qū)間均值,區(qū)間半徑

為了直觀對(duì)比修正前后的時(shí)域響應(yīng)區(qū)間，分別計(jì)算參數(shù)修正前后以及試驗(yàn)(真實(shí))的時(shí)域響應(yīng)區(qū)間，修正前后的自由端的時(shí)域響應(yīng)區(qū)間的修正值和試驗(yàn)值的對(duì)比，如圖9所示。圖9中的曲線代表的對(duì)象和“算例一”相同，從圖中可以看出兩種非線性形式下修正后的時(shí)域響應(yīng)區(qū)間的上下界和真實(shí)的響應(yīng)區(qū)間基本重合，證明兩種非線性形式下修正后的響應(yīng)區(qū)間都能反映真實(shí)的響應(yīng)區(qū)間。

圖9 修正前后時(shí)域響應(yīng)區(qū)間

4 結(jié) 論

本文基于區(qū)間攝動(dòng)原理，提出了一種求解不確定非線性結(jié)構(gòu)模型修正問題的方法。通過區(qū)間攝動(dòng)原理，將不確定性參數(shù)和時(shí)域動(dòng)響應(yīng)分別表示為區(qū)間均值和區(qū)間半徑的形式。通過建立不確定參數(shù)和響應(yīng)的區(qū)間均值、半徑映射關(guān)系將不確定性模型修正問題轉(zhuǎn)化為修正區(qū)間均值和區(qū)間半徑的確定性優(yōu)化問題。通過兩種不同非線性形式以及潛在線性結(jié)構(gòu)對(duì)比研究表明：

(1)立方非線性結(jié)構(gòu)和潛在的線性結(jié)構(gòu)的對(duì)比，兩種情況下修正后不確定參數(shù)的修正誤差均較小，說明此方法對(duì)非線性和線性結(jié)構(gòu)都能適用。

(2)初始參數(shù)區(qū)間對(duì)算法修正結(jié)果影響較低，改變初始的參數(shù)區(qū)間，修正的區(qū)間參數(shù)仍能快速收斂，且修正后的結(jié)果具有很高的精度。

(3)立方非線性和分段非線性兩種非線性結(jié)構(gòu)的不確定參數(shù)的區(qū)間均值和區(qū)間半徑都能快速的收斂到真實(shí)值附近，證明對(duì)于不同的非線性形式，該方法都能進(jìn)行高效的區(qū)間修正。

根據(jù)修正后的區(qū)間構(gòu)造參數(shù)樣本，代入有限元模型中計(jì)算對(duì)應(yīng)的響應(yīng)區(qū)間，修正之后的響應(yīng)區(qū)間與試驗(yàn)樣本(真實(shí)值)的響應(yīng)區(qū)間吻合良好，能夠準(zhǔn)確的反映試驗(yàn)樣本的區(qū)間特性。