柳振海,李 剛,劉 博,徐明強,呂 晴,蔣 上
1.(中國電力工程顧問集團有限公司,北京 100120;2.中國海洋大學(xué),山東 青島 266100 )
海上風電作為一種清潔能源,越來越受到世界各國的重視,大力發(fā)展海上風電,可優(yōu)化能源結(jié)構(gòu),促進經(jīng)濟、社會、環(huán)境協(xié)調(diào)發(fā)展。然而,海上風電設(shè)施長期服役于復(fù)雜的海洋環(huán)境中,一旦發(fā)生故障或損傷,不僅影響風電運行,也會造成巨大的經(jīng)濟損失,開展結(jié)構(gòu)健康監(jiān)測勢在必行。獲取準確的有限元模型是結(jié)構(gòu)健康監(jiān)測的重要內(nèi)容,但由于建模誤差、參數(shù)不確定性等問題,有限元模型不一定能準確反映結(jié)構(gòu)的真實情況,需要對有限元模型進行修正,使得修正后的模型能準確反映真實結(jié)構(gòu)的動力響應(yīng)特征。
有限元模型修正方法可以分為確定性模型修正和不確定性模型修正。前者將結(jié)構(gòu)參數(shù)及響應(yīng)假設(shè)為恒定不變的,修正結(jié)果為唯一確定解。但在實際工程設(shè)計建設(shè)中,結(jié)構(gòu)部件的實際加工組裝可能與設(shè)計值存在一定誤差,導(dǎo)致結(jié)構(gòu)參數(shù)的不確定性。因此確定性模型修正得到的最優(yōu)解只是不確定性模型修正解中的一個特例。隨機模型修正是對傳統(tǒng)模型修正方法的創(chuàng)新,并考慮到了不確定性因素的影響,利用統(tǒng)計概率分析量化不確定性因素,將不確定模型修正問題轉(zhuǎn)化為均值和標準差的修正問題,從而得到與實際結(jié)構(gòu)動力響應(yīng)特征一致的有限元模型[1-3]。
由于能夠代替復(fù)雜的有限元模型分析,代理模型被廣泛用于模型修正。許澤偉[4]等提出一種基于多項式混沌展開和KL 散度的隨機有限元模型修正方法,以三維桁架為例,對彈性模量和密度的均值和標準差進行修正。冷建成[5]等以某海洋平臺結(jié)構(gòu)為研究對象,提出了Kriging模型與多目標遺傳算法優(yōu)化相結(jié)合的動力學(xué)模型修正方法,實驗顯示能夠顯著改善Kriging 模型精度。孫永朋[6]等針對隨機模型修正精度和效率低的問題,提出一種基于Kriging 模型和小波包能量譜的隨機有限元模型修正方法。丁雅杰[7]等提出一種基于貝葉斯推理的非線性結(jié)構(gòu)模型修正方法,同時考慮激勵的隨機性,建立了復(fù)合隨機振動系統(tǒng)的動力可靠度分析方法。
本文針對小樣本數(shù)據(jù)的模型修正提出了一種基于KL 散度的分步型隨機模型修正方法,通過數(shù)值模擬驗證了此方法的有效性。首先采用靈敏度分析法選擇待修正參數(shù),通過正交試驗生成訓(xùn)練樣本,利用有限元模型計算樣本響應(yīng),從而構(gòu)造以待修正參數(shù)為輸入,以頻率響應(yīng)為輸出的Kriging 代理模型;利用蒙特卡洛模擬得到初始樣本點并計算其響應(yīng)值,以樣本點響應(yīng)與實測數(shù)據(jù)的KL 散度構(gòu)造目標函數(shù),通過多目標優(yōu)化算法對海上風電結(jié)構(gòu)進行隨機模型修正。
相較于其他響應(yīng)面模型(如多項式回歸、神經(jīng)網(wǎng)絡(luò)等),Kriging 模型在構(gòu)造過程中只需要較少的樣本點,計算成本較低,并且處理非線性數(shù)據(jù)效果優(yōu)異[8],能實現(xiàn)對有限元模型響應(yīng)較準確預(yù)測,更適合處理小樣本實測數(shù)據(jù)的模型修正。
Kriging 模型是一種通過已知樣本點信息估計未知試驗點信息的無偏估計模型,采用Kriging 代理模型的前提是假設(shè)所有數(shù)據(jù)之間都服從n維正態(tài)分布[9]。
已知樣本點x= [x1,x2, ...,xm],xi是n維向量,對應(yīng)響應(yīng)為y(x) = [y(x1),y(x2), ...,y(xm)],由于x服從正態(tài)分布,因此y服從均值為μ、協(xié)方差C的多維高斯分布。y(x)可以用式(1)表示:
式中:P代表基函數(shù)的個數(shù);f(x)=[f1(x),f2(x),…,fl(x),…,fp(x)]∈Rp×p為基函數(shù)矩陣;β=[β1,β2…,βl,…,βP]為基函數(shù)權(quán)重;fT(x)β為線性回歸模型,z(x)為服從高斯分布的隨機過程;不同變量之間的協(xié)方差矩陣C為:
式中:R為數(shù)據(jù)點之間的協(xié)方差矩陣;σ2為方差。
樣本點的似然概率為:
式中:n為樣本數(shù)。
對式(3)取對數(shù):
使得似然概率最大的超參數(shù)β,σ2的估計值和為:
對于建立的代理模型,還需要對其進行進一步的驗證。本文采用決定系數(shù)R2和相對均方誤差RMSE評價Kriging 代理模型精度,R2的取值范圍為0 ~1,R2值越大,表示預(yù)測效果越好,RMSE越趨近于0 表示預(yù)測效果越好,其表達式為
KL 散度可以用來度量兩個概率分布之間的差異[10],需要注意的是,KL 散度并不滿足對稱性,即DKL(P||Q) ≠DKL(Q||P)。假設(shè)實測頻率與樣本頻率分別為二維正態(tài)分布P和Q,此時DKL(P||Q)為前向KL 散度,主要表現(xiàn)為在任何P(x) > 0 的位置使得Q(x) > 0,即使得Q能夠最大程度上覆蓋P;DKL(Q||P)為反向KL散度,主要表現(xiàn)為在P(x)趨近于0 時,Q(x)也盡可能趨近于0。
在優(yōu)化問題中,需要通過迭代優(yōu)化樣本頻率Q使得兩個分布之間的KL 散度能夠盡可能小,從而實現(xiàn)用樣本頻率Q擬合實測頻率P,且實測頻率分布往往是多峰分布,前向KL 散度求解得到多個極小值對應(yīng)的待修正參數(shù)的平均值,不能滿足優(yōu)化要求,而采用反向KL 散度則可以確保至少達到一個局部最優(yōu)解,因此本文采用KL 散度計算實測頻率與樣本頻率分布之間的差異程度。
假設(shè)P、Q服從均值為μ0,μ1,協(xié)方差矩陣為Σ0,Σ1的高斯分布,則兩分布之間的KL散度為:
式中:tr表示矩陣跡;det 表示某矩陣的行列式;k為Σ0的維度。
在優(yōu)化算法中,目標函數(shù)提供了對不同解進行評估和比較的指標,決定了算法的搜索目標,因此選擇適當?shù)哪繕撕瘮?shù)對解決實際問題至關(guān)重要。本文以待修正參數(shù)的均值和方差為修正目標,以實測頻率與樣本點頻率的KL 散度為目標函數(shù),同時對參數(shù)的均值和方差進行修正。由于兩個分布沒有重疊時,KL 散度沒有意義,無法進行迭代操作,因此將目標函數(shù)分為兩個部分:首先通過樣本點頻率及實測頻率的均值、方差構(gòu)造多目標函數(shù)進行優(yōu)化計算,當兩個分布出現(xiàn)重疊時,改用KL 散度為目標函數(shù),具體目標函數(shù)構(gòu)造如下:
式中:nsim為樣本響應(yīng)數(shù)量;nobs為實測值數(shù)量。
在確定好目標函數(shù)后,需要選擇合適的優(yōu)化算法求待修正參數(shù)的最優(yōu)解。NSGA-Ⅱ[11]是目前最常見的多目標優(yōu)化算法之一,在原始算法的基礎(chǔ)上,引入了快速非支配排序、擁擠度等概念,減少了計算的復(fù)雜度,并提高了算法效率,其基本思想為:
1)隨機生成規(guī)模為N 的初始種群Pt,通過非支配排序、選擇、交叉、變異,生成子代種群Qt;
2)將初始種群Pt與子代種群Qt合并得到新的種群Rt;
3)進行快速非支配排序,并對每個非支配層的個體進行擁擠度計算,根據(jù)非支配關(guān)系及個體的擁擠度選擇合適的個體組成新的父代種群Pt+1;
4)采用遺傳算法生成新的子代種群Qt+1,并將Pt+1與Qt+1合并生成新的種群Rt;
5)依此類推,直到滿足程序結(jié)束的條件。
選取一個5 MW 海上風機作為研究對象,風機采用導(dǎo)管架基礎(chǔ),基礎(chǔ)深入海底,底部高程為-50.001 m,頂部高程為10.262 m,上部塔筒高度100 m,采用變截面,塔筒底部直徑6 m,壁厚0.035 1 m,頂部直徑3.87 m,壁厚0.024 7 m,頂部機箱和葉片簡化為偏心質(zhì)量點,共31.452 t。利用ANSYS 建立有限元模型,模型共166 個節(jié)點,221 個單元。對其進行模態(tài)分析,得到前6 階頻率為:0.490 16 Hz、0.490 28 Hz、1.661 4 Hz、1.661 8 Hz、3.774 2 Hz、3.774 3 Hz,前三階模態(tài)振型如圖1 所示。
圖1 海上風機前三階模態(tài)振型
考慮到海上風機長期受腐蝕、沖刷、溫度、濕度等因素影響,其中材料參數(shù)選擇彈性模量、密度為待選參數(shù),幾何參數(shù)選擇基礎(chǔ)腿柱外徑、斜撐外徑、腿柱壁厚、斜撐壁厚為待選參數(shù)。對待選參數(shù)進行靈敏度分析,選擇高靈敏度修正參數(shù):將各初選參數(shù)值增加2%,得到修改后的頻率值,計算各階頻率變化率,除以參數(shù)變化率即為各參數(shù)前兩階頻率的靈敏度S,如式(11)所示。
式中:f i為第i階初始頻率,Hz;為參數(shù)值增加后的第i階初始頻率,Hz;Pm為初始參數(shù)值;Pmc為增加后的參數(shù)值。
表1 為修正參數(shù)各階頻率靈敏度,由表1可知,在材料參數(shù)中,彈性模量與密度靈敏度相差不大,在幾何參數(shù)中,腿柱外徑靈敏度較高。由于采用材料參數(shù)可以有效地提高方差修正效果,因此最終彈性模量(E)、密度(Des)及腿柱外徑(OL)為待修正參數(shù)。假設(shè)待修正參數(shù)服從高斯分布,根據(jù)施工圖紙及工程經(jīng)驗確定其初始均值、方差見表2 所列。
表1 修正參數(shù)各階頻率靈敏度
表2 修正參數(shù)初始均值與方差
根據(jù)各參數(shù)的初始均值和方差選擇95%置信區(qū)間為修正參數(shù)取值范圍,即:彈性模量取值范圍為[205.941,206.059],密度取值范圍為[7 849.4,7 850.6],腿柱外徑取值區(qū)間為[1.191,1.309]。通過正交試驗設(shè)計選擇,取三因素五水平共25 組設(shè)計樣本參數(shù),見表3 所列,將樣本參數(shù)代入有限元模型中計算得到前兩階頻率,由此構(gòu)建Kriging 模型。
表3 正交試驗設(shè)計樣本
采用相對均方誤差RMSE與決定系數(shù)R2對建立的Kriging 代理模型進行有效性評價,計算結(jié)果見表4 所列,證明構(gòu)造的Kriging 代理模型精度較高,可以代替有限元模型。
表4 Kriging代理模型有效性評價
假設(shè)有限元模型待修正參數(shù)試驗值與初始值見表5 所列,對試驗值進行蒙特卡洛抽樣1 000 次得到樣本點,通過Kriging 模型計算得到樣本點前兩階頻率響應(yīng)作為實測數(shù)據(jù)。首先采用實測數(shù)據(jù)與每次迭代生成的隨機樣本頻率的均值和方差構(gòu)造多目標函數(shù),隨后利用NSGA-Ⅱ多目標優(yōu)化算法進行求解,得到結(jié)構(gòu)待修正參數(shù)的修正值,見表5 所列。
表5 參數(shù)修正前后誤差
表5 中標準差是實測數(shù)據(jù)與其平均數(shù)離差平方的算術(shù)平均數(shù)的平方根,反映實測數(shù)據(jù)的離散程度,誤差是試驗值與初始值之間的差異。根據(jù)修正參數(shù)最優(yōu)解,通過構(gòu)造的Kriging 模型進行計算,得到修正后有限元模型前五階頻率的均值和標準差,與初始值和試驗值進行比較,見表6、表7 所列,修正前后頻率均值的最大誤差由68.59%降低為0.1%,標準差的最大誤差由212.82%降低為17.29%,且前三階頻率標準差最大誤差為0.48%,修正結(jié)果表明本文方法有效。
表6 修正前后頻率均值
繪制修正前后模型與試驗?zāi)P皖l率置信橢圓如圖2 所示,可以看出,修正前實測值與初始值頻率置信橢圓大小與中心點位置有明顯差異,修正后實測值與修正值的頻率置信橢圓的大小、中心點、傾斜方向幾乎完全相同。圖3給出了結(jié)構(gòu)前兩階頻率修正前后的概率密度曲線,對比可知,修正前實測值曲線與初始值曲線的峰值有明顯差異,兩條曲線重合程度較低。修正后的實測值與修正值概率密度曲線形狀相似且?guī)缀跬耆呛?。進一步驗證了本文所提修正方法的有效性。
圖2 修正前后模型與試驗?zāi)P皖l率置信橢圓
圖3 修正前后頻率概率密度曲線
考慮到參數(shù)隨機不確定對結(jié)構(gòu)響應(yīng)分布的影響以及KL 散度在不確定性度量的優(yōu)勢,本文采用Kriging 模型代替有限元模型進行計算,并通過靈敏度分析法選擇合適的待修正參數(shù),提高了代理模型的精度和修正效率。同時,選擇能夠度量兩個分布之間差距的KL 散度構(gòu)造目標函數(shù),以海上風機頻率均值和標準差為修正目標進行隨機模型修正,修正后的有限元模型響應(yīng)與實際結(jié)構(gòu)響應(yīng)高度吻合。結(jié)果表明,基于Kriging 模型和KL 散度的隨機模型修正方法可以用于小樣本實測數(shù)據(jù)的海上風電結(jié)構(gòu)隨機模型修正,且修正效果顯著。開展風電實驗?zāi)P托拚龑⑹俏磥磉M行的工作。