基于響應(yīng)面全局優(yōu)化技術(shù)的蜂窩板材料性能參數(shù)修正

孫衛(wèi)青，程 偉

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191)

近年來，航天器的微振動問題已經(jīng)引起越來越多的關(guān)注。微振動是由各種擾振源(如動量輪，控制力矩陀螺和太陽能帆板作動器)所激勵引起，而后通過航天器的結(jié)構(gòu)傳遞到成像或指向裝置，導(dǎo)致精度降低。微振動的能量分布在1kHz左右的較寬的頻率范圍內(nèi)，可能會激發(fā)起航天器上某些結(jié)構(gòu)件的彈性體模態(tài)，而且這種寬頻的振動又無法通過姿態(tài)調(diào)整系統(tǒng)或軌道控制系統(tǒng)來進行抑制。

對于微振動的控制有兩種途徑。一方面，國內(nèi)外有很多學(xué)者在研究通過設(shè)計無源或有源隔振器來最小化擾振源作用到結(jié)構(gòu)上的載荷。而另一方面，可以通過優(yōu)化航天器結(jié)構(gòu)，并采用具有高剛度/質(zhì)量比的材料(例如目前廣泛使用的蜂窩夾層板)，從而降低傳遞到最終目標點的振動。

對于蜂窩夾層板結(jié)構(gòu)的動力學(xué)特性研究，早在20世紀80年代，基于Mindlin板殼理論，曾經(jīng)開發(fā)出了一些夾層板彎曲振動的理論方法[1-4]。而自20世紀90年代以來，有限元(FE)方法開始廣泛應(yīng)用于結(jié)構(gòu)動力學(xué)仿真[5-13]。無論采用何種夾層板結(jié)構(gòu)的動力學(xué)分析方法，要想獲取準確的分析結(jié)果，都首先要保證表層及特別是芯層結(jié)構(gòu)的材料參數(shù)是精確的。在上述的一些研究工作中，芯層的材料參數(shù)往往都是通過實驗測定的。然而，對于一些軟芯材料(如蜂窩芯)，因為在實驗中往往需要通過特別的夾緊裝置對其進行固定，因而想要準確地獲取其材料參數(shù)尤其是面外剪切模量的難度較大。而芯層的面外剪切模量又恰恰對于整體夾層結(jié)構(gòu)的動態(tài)特性而言是最為重要的。由于構(gòu)建詳細的芯層有限元模型計算成本較高，因此工程上最有效的方法是基于理論計算構(gòu)建芯層甚至整個夾層板的等效模型，但是考慮到生產(chǎn)加工過程中的工藝誤差，必須對等效模型進行驗證和校核。Jiang等[14]曾試圖通過模態(tài)實驗數(shù)據(jù)對鋁蜂窩板的蜂窩芯的面外剪切模量進行修正。Debruyne等[15]通過對模態(tài)實驗數(shù)據(jù)及有限元分析結(jié)果進行隨機性分析，研究了表層蒙皮楊氏模量和蜂窩芯面外剪切模量的變化特性。孔憲仁等[16]曾試圖將芯層等效為各向同性材料，建立了材料參數(shù)對應(yīng)蜂窩板固有頻率的響應(yīng)面模型，但并沒有通過實驗進行有限元模型的有效性驗證?？傮w而言，目前針對這方面的研究工作并不多，而且對于芯層材料等效參數(shù)辨識的優(yōu)化算法研究也不甚完整。

對于航天器的微振動問題，由于其本身具有激勵能量小，振動響應(yīng)級別在微米級別的特點，因此對于有限元模型的精度有更高的要求。

此外，航天器蜂窩板結(jié)構(gòu)中通常有用于裝配的預(yù)埋件，對于預(yù)埋件的處理目前鮮有文獻提及。本工作建立了典型航天器蜂窩板帶預(yù)埋件結(jié)構(gòu)的動力學(xué)模型，并考慮了膠層附加質(zhì)量的影響，嘗試使用基于響應(yīng)面模型的快速全局優(yōu)化技術(shù)對蜂窩板動力學(xué)有限元模型進行修正，以滿足微振動分析的高精度要求。

1 方法及原理

1.1 雙層壁蜂窩芯的等效模型

給定芯層材料的楊氏模量Es，剪切模量Gs，密度ρ及泊松比νs，可以根據(jù)以下理論計算雙厚度壁六邊形蜂窩芯(如圖1所示)的等效材料特性參數(shù)。

圖1 雙層壁厚六面體蜂窩芯胞元示意圖Fig.1 Sketch of a cell in a hexagonal honeycomb core

首先通過施加單位載荷和單位位移的方法推導(dǎo)出兩個方向的面外剪切剛度值G13和G23[17]

(1)

(2)

根據(jù)胞元中材料所占的百分比，可以計算得到芯層材料的密度

(3)

考慮蜂窩芯的法向剛度與其密度成正比，得到

(4)

通過在末端施加力矩，然后計算胞壁的彈性變形，獲得面內(nèi)模量和泊松比[18-19]

(5)

(6)

(7)

(8)

(9)

基于上述方程的推導(dǎo)過程，進一步綜合考慮了軸向和剪切變形，推導(dǎo)出了精確的面內(nèi)泊松比[20]

(10)

對于剪切模量G23，可以通過對蜂窩芯單胞進行有限元分析并對比不同高寬比的剪切模量，得到G23的最小二乘回歸近似確定值[21]

(11)

面外泊松比ν31和ν32等于芯層材料的泊松比，根據(jù)互易性關(guān)系計算得到ν13和ν23[22]

ν31=ν32=νs

(12)

(13)

(14)

1.2 響應(yīng)面模型

響應(yīng)面模型(response surface model, RSM)目前已經(jīng)在很多領(lǐng)域得以廣泛應(yīng)用[23]。王開山等[24]曾嘗試用二階多項式構(gòu)建了印制電路板的動力學(xué)響應(yīng)面模型。鮑諾等[25]針對某飛機模型的前10階模態(tài)頻率建立了四階多項式響應(yīng)面模型。

要得到一個準確的響應(yīng)面模型，需要預(yù)先考慮兩個方面：數(shù)值分析樣本的實驗設(shè)計方法及構(gòu)建響應(yīng)面模型的函數(shù)算法。

1.2.1 實驗設(shè)計方法

實驗設(shè)計的抽樣方法大致可分為兩類：正交設(shè)計和隨機設(shè)計。正交設(shè)計的優(yōu)點是所有輸入?yún)?shù)在統(tǒng)計學(xué)上都是獨立的，因此在物理實驗中隨機誤差最小。然而，正交設(shè)計總是集中在邊界上[26]，實驗次數(shù)將隨著輸入?yún)⒘吭O(shè)計值的數(shù)量成指數(shù)倍數(shù)增長，對于高階響應(yīng)面模型而言將會大幅提高計算成本。因此，正交設(shè)計一般主要用于設(shè)計參量篩選，從大量設(shè)計變量中識別出最重要的因素[27]。

對于確定性計算分析(包括有限元分析)不存在隨機誤差，而且隨機采樣可以最小化設(shè)計區(qū)域上的累計均方差[28]。目前最常用的隨機抽樣方法是拉丁超立方體設(shè)計方法(Latin hypercube design，LHD)。

圖2是中心復(fù)合實驗設(shè)計(central composite ins-cribed design，CCID)和一種LHD對于兩個輸入?yún)⒘?X1，X2)的抽樣對比。

圖2 中心復(fù)合實驗設(shè)計(a)和一種拉丁超立方體實驗設(shè)計(b)對于兩個輸入?yún)⒘康某闃訉Ρ菷ig.2 Sample points for two input factors of CCID(a) and a possible LHD(b)

1.2.2 響應(yīng)面模型構(gòu)建算法

最常用的響應(yīng)面模型構(gòu)建算法包括多項式回歸算法、Kriging算法和徑向基函數(shù)(RBF)等算法。Jin等[29]通過比較不同模型對于具有不同問題尺度，不同程度非線性和噪聲的問題所獲得的結(jié)果，認為在大多數(shù)情況下，RBF是在準確度和魯棒性方面最可靠的方法，特別是對于高階非線性問題而言。

RBF模型的基函數(shù)以采樣數(shù)據(jù)點和待預(yù)測點之間的歐幾里德距離為變量，這樣可以保證RBF模型在所有采樣點位置的函數(shù)值等于采樣數(shù)據(jù)的輸出響應(yīng)。RBF模型的一般形式是

(15)

(16)

其中βi是加權(quán)因子，ri表示輸入向量和第i次實驗設(shè)計點之間的歐幾里德距離，φ(r)是歐幾里德距離函數(shù)，其中最為普遍的函數(shù)形式包括以下4種：

(1)線性函數(shù)：φ(r)=rφ(r)=r

(2)立方函數(shù)：φ(r)=r3

(3)薄板樣條：φ(r)=r2lg(r)

基于實驗樣本，方程(15)可以以矩陣形式寫成

Aβ=y

(17)

其中，

(18)

β=(β1,…,βn)T

(19)

y=(y(1),…,y(n))T

(20)

通過線性方程組求解來，可以確定得到各加權(quán)因子βi(i= 1，…，n)

β=A-1y

(21)

構(gòu)造RBF模型完成后，就可以模擬輸入?yún)?shù)空間內(nèi)任意點x的輸出

(22)

1.3 基于響應(yīng)面模型的優(yōu)化算法

基于響應(yīng)面模型的優(yōu)化算法是一種基于響應(yīng)面模型進行自適應(yīng)采樣的全局優(yōu)化算法。Jones等[30]最早提出了基于Kriging響應(yīng)面模型進行自適應(yīng)LHD采樣進行全局優(yōu)化的方法，其原理如圖3所示。首先使用基于LHD采樣的少數(shù)幾個樣本來初步探索設(shè)計空間，然后建立初步的響應(yīng)面模型，根據(jù)得到的響應(yīng)面模型分析可能的最優(yōu)區(qū)域，在該區(qū)域中逐步添加新的樣本，并不斷更新響應(yīng)面模型。重復(fù)該過程，直到新的采樣點非常接近現(xiàn)有采樣點之一。優(yōu)化終止條件通?？梢远x為新采樣點與現(xiàn)有采樣點間的距離容差，即當新采樣點與現(xiàn)有采樣點的距離相對現(xiàn)有采樣點空間坐標的百分比小于該容差限值時，優(yōu)化終止。

評估添加新樣本的可能最優(yōu)區(qū)域主要基于兩方面的考慮：(1)響應(yīng)面模型是否能夠預(yù)測到更好的設(shè)計；(2)在該區(qū)域內(nèi)的樣本數(shù)是否太少以至于響應(yīng)面模型在該區(qū)域的準確性無法保證。因為這種優(yōu)化算法僅基于自適應(yīng)采樣，相比較遺傳基因等全局優(yōu)化算法的迭代過程，可以充分保證其效率。

圖3 基于響應(yīng)面模型的優(yōu)化算法流程圖Fig.3 Efficient global optimization work flow for response surface modal

2 實驗及有限元建模

圖4所示的蜂窩夾心板結(jié)構(gòu)中，面板和芯層材料均為鋁合金(材料特性如表1所示)。根據(jù)式(1)～(14)計算得到蜂窩芯的等效材料性能(如表2所示)。

圖4 蜂窩板尺寸示意圖Fig.4 Schematic diagram of honeycomb plate

Young’s modulus,E/GPaDensity, ρ/(kg·m-3)Poissonratio,ν69 2730 0.33

如圖5所示，通過錘擊法模態(tài)實驗測量夾層板的動態(tài)響應(yīng)特性。被測蜂窩板用軟橡膠帶懸掛以模擬自由邊界條件。總共有56個敲擊點(7×8)和1個加速度響應(yīng)點，帶寬設(shè)置為2kHz，頻率分辨率為0.5Hz。相關(guān)測試設(shè)備列于表3中。

表2 蜂窩芯等效材料特性參數(shù)Table 2 Equivalent material properties of the honeycomb core

圖5 蜂窩板動態(tài)響應(yīng)實驗示意圖Fig.5 Test set-up of dynamic response experiment for the honeycomb plate

表3 模態(tài)測試設(shè)備列表Table 3 Measurement devices employed for modal testing

使用Nastran以殼-實體-殼(shell-volume-shell，SVS)的方式進行蜂窩板建模，其中兩個面板使用殼單元，芯層使用實體單元。整個模型包含5376個實體單元和2層共5376個殼單元。

對于用于螺栓連接的鋁材質(zhì)預(yù)埋件，在模型中處理為分布式集中質(zhì)量添加到對應(yīng)的芯層實體節(jié)點上。在整個蜂窩板中有兩種預(yù)埋件，兩個垂直邊上有12個大預(yù)埋件，每個質(zhì)量為0.0405kg，兩個水平邊上有14個小預(yù)埋件，每個質(zhì)量為0.0183kg。膠層的附加質(zhì)量影響也作為分布式集中質(zhì)量附加到兩層面板的殼單元節(jié)點上。膠層的質(zhì)量由整板實際質(zhì)量減去鋁蜂窩芯、兩個鋁面板及預(yù)埋件質(zhì)量而計算得到，總質(zhì)量為0.60kg。

表4為通過有限元分析獲得的1kHz帶寬內(nèi)前6個模態(tài)頻率與實驗結(jié)果的對比。各階模態(tài)的模態(tài)置信準則(MAC)值均大于0.9，最大頻率誤差不超過3%，平均誤差為1.24%。說明有限元分析結(jié)果與實驗結(jié)果吻合較好。但對于第1階、第2階以及第5階模態(tài)，誤差在1%～3%之間，可期待通過模型修正進一步提高其精度。

表4 初始有限元模型分析結(jié)果與實驗結(jié)果模態(tài)頻率對比Table 4 Modal frequencies from the FE analysis and experimental results

3 有限元模型修正

3.1 識別用于優(yōu)化的輸入?yún)?shù)

在進行有限元模型修正之前，首先需要確定待優(yōu)化參數(shù)。對于蜂窩板蒙皮的材料參數(shù)、預(yù)埋件及膠層的質(zhì)量均可認為是確定的?？紤]到蜂窩板加工工藝的誤差，應(yīng)將蜂窩芯理論等效模型的材料屬性作為待優(yōu)化參數(shù)。

蜂窩芯屬于正交異性材料，對應(yīng)Nastran實體單元CHEXS，需要將蜂窩芯的正交各向異性材料參數(shù)根據(jù)胡克定律轉(zhuǎn)移為MAT9的各向異性材料參數(shù)矩陣。即將6個材料參數(shù)(E11，E22，E33，G12，G23，G31)和ρ轉(zhuǎn)化到6×6的MAT9材料參數(shù)矩陣D中，矩陣D中8個有效元素與6個材料參數(shù)的對應(yīng)關(guān)系為：d11=(1-ν23ν32)/(E22E33Δ)，d12=d21=(ν12+ν13ν32)/(E11E33Δ)，d13=d31=(ν31+ν21ν32)/(E22E33Δ)，d22=(1-ν13ν31)/(E11E33Δ)，d33=(1-ν12ν21)/(E11E22Δ)，d44=G12，d55=G23，d66=G31，其中Δ=(1-ν12ν21-ν23ν32-ν31ν13-2ν12ν23ν31)/(E11E22E33)，其他元素均為零。

鑒于輸入?yún)⒘肯鄬^多，為提高優(yōu)化效率，可以首先確定出仿真模型對哪些輸入?yún)?shù)最為敏感，哪些是最為重要的輸入?yún)?shù)，以縮減優(yōu)化變量的個數(shù)。

利用1.2節(jié)中提到的正交實驗設(shè)計中的中心復(fù)合實驗設(shè)計(CCID)進行計算，然后根據(jù)計算產(chǎn)生的樣本進行Pearson相關(guān)分析，可以確定出最重要的輸入?yún)?shù)。

輸出結(jié)果(y)和第k個輸入?yún)?shù)(xk)之間的Pearson相關(guān)系數(shù)rk定義為：

(23)

Pearson相關(guān)系數(shù)取值在-1和+1之間。接近+1或-1意味著該輸入?yún)?shù)對分析結(jié)果存在直接的正向或反向影響，而接近0表示這個輸入?yún)?shù)與分析結(jié)果之間相對獨立。相關(guān)性分析的樣本能夠覆蓋整個設(shè)計空間，因此它可以反映每個輸入因子對結(jié)果的總體影響。

表5列出了各輸入?yún)⒘繉?yīng)各階固有頻率的相關(guān)系數(shù)。圖6對比了各參量對六階模態(tài)頻率的平均相關(guān)系數(shù)，可以確定出d55(G23)，d66(G31)和ρ是3個最重要的參數(shù)。

3.2 基于響應(yīng)面模型的優(yōu)化

根據(jù)相關(guān)性分析的結(jié)果，選擇蜂窩芯的3個等效材料參數(shù)(G23，G31和ρ)作為進一步優(yōu)化的輸入?yún)?shù)，將有限元分析和實驗結(jié)果的前六階模態(tài)頻率的平均誤差作為最小化目標函數(shù)：

表5 蜂窩芯材料參數(shù)對應(yīng)各階模態(tài)的Pearson相關(guān)系數(shù)Table 5 Correlation factors of the honeycomb core material parameters for each modal frequency

圖6 蜂窩芯各材料參數(shù)與六階模態(tài)頻率的平均相關(guān)系數(shù)Fig.6 Averaged correlation factors for the six modal frequencies

(24)

其中λi是各階模態(tài)特征頻率，下標i是模態(tài)階數(shù)。

使用徑向基函數(shù)響應(yīng)面模型進行自適應(yīng)拉丁超立方采樣，進行全局優(yōu)化。徑向基函數(shù)采用立方函數(shù)，每次迭代的采樣點數(shù)為10。優(yōu)化終止條件中新采樣點與現(xiàn)有采樣點的距離容差設(shè)為0.1%。圖7為整個迭代過程中的采樣點分布，其中不同顏色代表密度的大小(范圍為20～35kg/m3)，圓的直徑尺寸代表頻率誤差的大小(范圍為0.88%～22.24%)，從圖7中可以看到當G23= 130MPa，G31= 91MPa附近優(yōu)化過程得以收斂。整個優(yōu)化過程僅用了22次計算即到了輸入?yún)?shù)最優(yōu)解(見表6)。

圖7 基于徑向基函數(shù)響應(yīng)面模型的全局優(yōu)化采樣點分布Fig.7 Samples distribution during global optimization based on RBF response surface model

表6 蜂窩芯等效材料特性參數(shù)優(yōu)化結(jié)果Table 6 Optimization results of equivalent material property of the honeycomb core

表7為修正后有限元模型的模態(tài)頻率與實驗結(jié)果的對比，圖8為初始有限元模型和修正后有限元模型的模態(tài)頻率精度對比。從中可以看出除第3階，其他各階分析精度相比較初始有限元模型均得以提高，平均頻率誤差從1.24%降低為0.88%。而第3階模態(tài)誤差盡管略有所放大，但精度仍可控制在0.15%，相比較其他階，尤其是第4階和第5階模態(tài)精度的改善，是可以接受的。

表7 修正后有限元模型結(jié)果與實驗結(jié)果的模態(tài)頻率對比Table 7 Modal frequencies from the updated FE analysis and experimental results

圖9為當密度為23.58 kg/m3時，G23和G31相對模態(tài)頻率平均誤差的響應(yīng)面模型。圖9中可以看出，當2個面外剪切模量G23和G31取值分別超過100MPa和70MPa后，平均頻率誤差精度逐漸穩(wěn)定，逐漸逼近最優(yōu)點。

圖8 初始和修正后有限元模型的模態(tài)頻率精度對比Fig.8 Modal frequencies error comparison between initiative and updated FE model

圖9 G23和G31相對6階平均模態(tài)頻率誤差的響應(yīng)面 模型(ρ=23.58kg/m3)Fig.9 Response surface model ofG23andG31for average error of first six modal frequencies (ρ=23.58kg/m3)

4 結(jié)論

(1)基于蜂窩夾層板的蜂窩芯等效參數(shù)模型，建立了考慮安裝預(yù)埋件及膠層附加質(zhì)量的鋁蜂窩夾層板有限元模型。

(2)通過正交數(shù)值實驗設(shè)計(中心復(fù)合實驗設(shè)計)從蜂窩芯層的9個等效材料參數(shù)中篩選出3個重要參數(shù)作為有限元模型修正的輸入?yún)⒘俊?/p>

(3)利用基于響應(yīng)面模型的全局優(yōu)化技術(shù)對蜂窩芯的兩個面外剪切模量及密度進行了參數(shù)優(yōu)化，修正后的有限元模型的前六階模態(tài)的平均頻率誤差達到0.88%，最大誤差不超過2.5%，精度可以滿足航天器微振動分析的要求。