王春水，李大磊，張二亮

（鄭州大學(xué) 機械工程學(xué)院，鄭州 450001）

0 引言

實際結(jié)構(gòu)材料屬性的復(fù)雜性以及機械加工過程中的各種隨機誤差等因素會造成結(jié)構(gòu)參數(shù)的不確定性，使得材料相同、按相同方法加工出的名義尺寸一致的一組結(jié)構(gòu)件之間存在不確定性，其動力學(xué)特性相似但不盡形同，這就需要建立一個有效的模型來表征結(jié)構(gòu)的不確定性。樊征兵等[1]采用了隨機因子法、區(qū)間因子法和攝動法對彈性模量、密度和截面積參數(shù)存在不確定的桁架結(jié)構(gòu)進行模態(tài)分析，其中攝動方法適用參數(shù)較小的變動，對于一般情況蒙特卡羅模擬有很好的適用性，但其計算量很大，而響應(yīng)面模型可很好地提高效率簡化計算，F(xiàn)ang S E等[2]運用二次多項式響應(yīng)面模型和蒙特卡羅模擬的方法來量化隨機模型彈性模量、密度、慣性矩等參數(shù)的不確定性。而實際上引起結(jié)構(gòu)不確定性的參數(shù)有很多，這些參數(shù)在結(jié)構(gòu)空間上隨機分布，為更好地反映實際情況，把結(jié)構(gòu)參數(shù)看成是隨空間位置而變化的隨機變量的隨機場理論被引入工程結(jié)構(gòu)數(shù)據(jù)分析以描述材料參數(shù)的空間相關(guān)性。白長青等[3]把轉(zhuǎn)軸的彈性模量、剪切模量、密度等參數(shù)表示為一維隨機場函數(shù)，研究了不確定性轉(zhuǎn)子系統(tǒng)的隨機有限元建模及響應(yīng)分析；Feng X等[4]運用均勻隨機場表征各有限元單元彈性模量的不確定性，研究了簡支梁和層流板的動力學(xué)特性。

目前關(guān)于結(jié)構(gòu)不確定性的表征及參數(shù)修正主要是針對某些結(jié)構(gòu)或材料參數(shù)，如彈性模量、慣性矩、厚度等的參數(shù)模型進行的，對其均值及方差做不確定性分析，其模型本身的復(fù)雜度不高，雖然可以采用各種有效的方法進行隨機有限元分析，但并不能很精確地反映實際結(jié)構(gòu)。

針對結(jié)構(gòu)彈性模量的不確定性，提出一種基于隨機場理論的超參數(shù)模型來表征一組結(jié)構(gòu)屬性參數(shù)的不確定性，并用參數(shù)估計的方法對模型參數(shù)進行量化。為表征和量化這種不確定性，首先提出一個基于高斯核函數(shù)的隨機場模型來表征參數(shù)的不確定性，用K-L展開來模擬該隨機場超參數(shù)模型，對模型進行頻率分析，用多元核密度估計估計出輸出響應(yīng)的概率分布密度函數(shù)，把試驗和模擬結(jié)果進行對比運用最大似然估計來量化模型的參數(shù)。最后通過對一個算例進行數(shù)值試驗驗證了提出方法的有效性；同時考慮了結(jié)構(gòu)尺寸誤差對參數(shù)量化的影響。

1 不確定性的隨機場表征

機械加工、熱處理及材料自身的不確定性等誤差都會導(dǎo)致結(jié)構(gòu)的不確定性，這些誤差一般都很小且相互獨立，根據(jù)中心極限定理相互獨立的大量微小隨機變量，其分布服從高斯分布，高斯分布假設(shè)便于計算使問題可解。因此可采用高斯隨機場模型來描述一組結(jié)構(gòu)彈性模量的不確定性。

1.1 高斯隨機場模型

高斯隨機場有以下兩個特點[5]：1）其數(shù)學(xué)期望和方差為一個與位置坐標無關(guān)的常數(shù)，即為一隨機數(shù)，xi表示空間一點；2）其自協(xié)方差函數(shù)只與隨機場上兩點相對位置距離有關(guān)而與兩點的絕對位置坐標無關(guān)，即為兩點之間的相對距離，為自相關(guān)函數(shù)，相關(guān)距離L是其重要的參數(shù)，表示在相關(guān)距離內(nèi)的兩點的參數(shù)具有明顯的關(guān)聯(lián)性。建立隨機場的關(guān)鍵是構(gòu)建協(xié)方差矩陣。

在有限元框架下，要將連續(xù)隨機場離散為隨機變量向量。對于高斯隨機場，常用的離散方法有：譜表示法、Karhunen -Loeve展開法等。這里選用K-L展開法。

1.2 Karhunen-Loeve展開

K-L展開的實質(zhì)是將隨機場分解為一系列不相關(guān)的隨機變量和確定系數(shù)（特征函數(shù)、特征值）。其具有以下優(yōu)越性：對于任意類型隨機場都是均方收斂的；相較于其他方法，展開相同有限項時具有最小的均方誤差[6]。隨機場(,)EXω的K-L展開式為[7]：

式中：X= (x,y)為空間點坐標；為隨機場的期望： ξi(X)為零均值互不相關(guān)的高斯隨機序列；ω為隨機事件； λi和 φi(X)分別為隨機場協(xié)方差矩陣CE(x,y)的特征值和特征函數(shù)。

可見，隨機場K-L展開的關(guān)鍵是獲得協(xié)方差矩陣的特征值和特征函數(shù)。所研究的隨機場的協(xié)方差矩陣是定義在規(guī)則幾何空間域上的，其表達式易知，可較容易地得到其特征值和特征向量，從而可以較容易地實現(xiàn)隨機場的K-L展開。

2 結(jié)構(gòu)不確定性量化

對于所提出的高斯隨機場模型，表征結(jié)構(gòu)不確定性的參數(shù)是未知的，需要對其進行量化。將試驗和模擬的輸出響應(yīng)的分布特征進行對比，對其結(jié)果做最大似然估計來量化模型的參數(shù)。最大似然估計是一種在給定統(tǒng)計特性下估計參數(shù)的方法，而模型的輸出響應(yīng)的分布特征是未知的，這就需要對輸出響應(yīng)的概率分布做非參數(shù)估計，核密度估計是常用的一種非參數(shù)估計方法。

2.1 核密度估計

設(shè)S1,S2,...,Sn是獨立同分布隨機變量，其分布密度函數(shù)的多元核密度估計可表示為[8]：

2.2 最大似然估計

對結(jié)構(gòu)的輸出響應(yīng)做最大似然估計，估計出滿足概率分布密度的參數(shù)，把可能性最大的那個參數(shù)作為真實的估計量。對數(shù)似然函數(shù)更易處理，故對于如式(2)所示的核密度估計函數(shù)，其對數(shù)似然函數(shù)可以表示為：

其中，Qexp為試驗輸出響應(yīng)；Qmod為模擬輸出響應(yīng)；num為試驗點數(shù)；選擇多元高斯分布的核函數(shù)假設(shè)[8]，取窗寬為，核函數(shù)為

3 算例驗證

用MATLAB軟件進行數(shù)值試驗，試驗?zāi)Ｐ蜑閼冶哿海鋮?shù)為：長1m，橫截面積0.06m×0.06m，質(zhì)量密度為7800Kg/m3，有限元建模時將梁劃分為50個單元，如圖1所示。

圖1 懸臂梁模型

每個單元的彈性模量可以通過式(1)按梁的有限元單元展開得到。其中n=50，均值，協(xié)方差函數(shù)取為。則彈性模量E的不確定性隨機場就由標準差和相關(guān)距離L進行表征，超參數(shù)也即是試驗設(shè)計中的輸入?yún)?shù)，通過數(shù)值試驗進行量化。輸出響應(yīng)為每根梁的前四階固有頻率，這些梁的輸出響應(yīng)組成固有頻率矩陣S。取輸入?yún)?shù)為的500根梁的固有頻率矩陣Sexp為試驗數(shù)據(jù)，取輸入?yún)?shù)為σ∈{1.0,1.2,1.4,…,3}，L∈{10,20,30,…,100}的10000根梁的固有頻率矩陣Smod為模擬數(shù)據(jù)，對Sexp和Smod進行核密度估計和最大似然估計，驗證超參數(shù)的可辨識性和隨機場模型的有效性。

3.1 實驗數(shù)據(jù)分析

圖2 梁單元彈性模量的隨機分布

3.2 一維核密度估計

圖3 一維核密度估計

可見實驗數(shù)據(jù)的概率密度分布規(guī)律與L0=50和σ=2的實驗數(shù)據(jù)的概率密度分布規(guī)律最為接近。

3.3 多維核密度估計及最大似然估計

為了更加充分地說明問題，要對Sexp和Smod進行四維核密度估計并在此基礎(chǔ)上做關(guān)于參數(shù)的最大似然估計。由于一些點處的概率密度函數(shù)值很小，為避免連乘時似然函數(shù)為0，故先對概率密度函數(shù)值取對數(shù)再連加，同時為了使結(jié)果凸顯便于觀察，對對數(shù)似然函數(shù)值做指數(shù)處理。最終可得到似然函數(shù)的三維圖如圖4所示。

圖4 似然函數(shù)數(shù)值

4 結(jié)論

提出的基于隨機場理論的超參數(shù)模型能夠表征結(jié)構(gòu)的不確定性，反映結(jié)構(gòu)參數(shù)的不確定性與空間關(guān)聯(lián)性，可以描述實際結(jié)構(gòu)參數(shù)的隨機變化特性。算例表明，模型超參數(shù)可以通過結(jié)構(gòu)的輸出響應(yīng)進行量化，得到結(jié)構(gòu)超參數(shù)的標準差和自相關(guān)距離，具有一定的應(yīng)用性。提出的方法也適用于其他參數(shù)，也可對其他材料參數(shù)或結(jié)構(gòu)參數(shù)建立隨機場超參數(shù)模型并進行量化。

[1] 樊征兵,何歡,陳國平.不確定結(jié)構(gòu)的模態(tài)和動響應(yīng)分析[A].第十屆全國振動理論及應(yīng)用學(xué)術(shù)會議論文集(2011)上冊[C].2011.

[2] FANG S E, REN W X,PERERA R. A stochastic model updating method for parameter variability quantification based on response surface models and Monte Carlo simulation[J].Mechanical Systems and Signal Processing,2012(33)：83-96.

[3] 白長青,張紅艷.不確定性轉(zhuǎn)子系統(tǒng)的隨機有限元建模及響應(yīng)分析[J].動力學(xué)與控制學(xué)報,2012,10(3)：283-288.

[4] FENG X, LU Z X, YANG Z Y, et al. Analysis on the variances of material and structural properties based on random field theory[J]. Probabilistic Engineering Mechanics.2011(26)：222-230.

[5] 張繼周,繆林昌.基于隨機場理論的地基概率沉降分析[J].巖土工程學(xué)報,2010,32(7)：1059-1064.

[6] 李少龍,楊金忠,張家發(fā).KL展開在滲流場隨機分析中的初步應(yīng)用[J].長江科學(xué)院院報,2009,26(10)：39-43.

[7] BETZ W, PAPAIOANNOU I, STRAUD D. Numerical methods for the discretization of random fields by means of the Karhunen-Loève expansion[J].Computer Methods Applied Mechanics and Engineering.2014(271)：109-129.

[8] SOIZE C. A computational inverse method for identification of non-Gaussian random fields using the Bayesian approach in very high dimension[J].Computer Methods Applied Mechanics and Engineering,2011(07).