基于響應面全局優化技術的蜂窩板材料性能參數修正

孫衛青，程 偉

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

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

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

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

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

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

1 方法及原理

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

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

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

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

(1)

(2)

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

(3)

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

(4)

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

(5)

(6)

(7)

(8)

(9)

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

(10)

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

(11)

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

ν31=ν32=νs

(12)

(13)

(14)

1.2 響應面模型

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

要得到一個準確的響應面模型，需要預先考慮兩個方面：數值分析樣本的實驗設計方法及構建響應面模型的函數算法。

1.2.1 實驗設計方法

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

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

圖2是中心復合實驗設計(central composite ins-cribed design，CCID)和一種LHD對于兩個輸入參量(X1，X2)的抽樣對比。

圖2 中心復合實驗設計(a)和一種拉丁超立方體實驗設計(b)對于兩個輸入參量的抽樣對比Fig.2 Sample points for two input factors of CCID(a) and a possible LHD(b)

1.2.2 響應面模型構建算法

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

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

(15)

(16)

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

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

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

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

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

Aβ=y

(17)

其中，

(18)

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

(19)

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

(20)

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

β=A-1y

(21)

構造RBF模型完成后，就可以模擬輸入參數空間內任意點x的輸出

(22)

1.3 基于響應面模型的優化算法

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

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

圖3 基于響應面模型的優化算法流程圖Fig.3 Efficient global optimization work flow for response surface modal

2 實驗及有限元建模

圖4所示的蜂窩夾心板結構中，面板和芯層材料均為鋁合金(材料特性如表1所示)。根據式(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所示，通過錘擊法模態實驗測量夾層板的動態響應特性。被測蜂窩板用軟橡膠帶懸掛以模擬自由邊界條件?？偣灿?6個敲擊點(7×8)和1個加速度響應點，帶寬設置為2kHz，頻率分辨率為0.5Hz。相關測試設備列于表3中。

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

圖5 蜂窩板動態響應實驗示意圖Fig.5 Test set-up of dynamic response experiment for the honeycomb plate

表3 模態測試設備列表Table 3 Measurement devices employed for modal testing

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

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

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

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

3 有限元模型修正

3.1 識別用于優化的輸入參數

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

蜂窩芯屬于正交異性材料，對應Nastran實體單元CHEXS，需要將蜂窩芯的正交各向異性材料參數根據胡克定律轉移為MAT9的各向異性材料參數矩陣。即將6個材料參數(E11，E22，E33，G12，G23，G31)和ρ轉化到6×6的MAT9材料參數矩陣D中，矩陣D中8個有效元素與6個材料參數的對應關系為：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)，其他元素均為零。

鑒于輸入參量相對較多，為提高優化效率，可以首先確定出仿真模型對哪些輸入參數最為敏感，哪些是最為重要的輸入參數，以縮減優化變量的個數。

利用1.2節中提到的正交實驗設計中的中心復合實驗設計(CCID)進行計算，然后根據計算產生的樣本進行Pearson相關分析，可以確定出最重要的輸入參數。

輸出結果(y)和第k個輸入參數(xk)之間的Pearson相關系數rk定義為：

(23)

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

表5列出了各輸入參量對應各階固有頻率的相關系數。圖6對比了各參量對六階模態頻率的平均相關系數，可以確定出d55(G23)，d66(G31)和ρ是3個最重要的參數。

3.2 基于響應面模型的優化

根據相關性分析的結果，選擇蜂窩芯的3個等效材料參數(G23，G31和ρ)作為進一步優化的輸入參數，將有限元分析和實驗結果的前六階模態頻率的平均誤差作為最小化目標函數：

表5 蜂窩芯材料參數對應各階模態的Pearson相關系數Table 5 Correlation factors of the honeycomb core material parameters for each modal frequency

圖6 蜂窩芯各材料參數與六階模態頻率的平均相關系數Fig.6 Averaged correlation factors for the six modal frequencies

(24)

其中λi是各階模態特征頻率，下標i是模態階數。

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

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

表6 蜂窩芯等效材料特性參數優化結果Table 6 Optimization results of equivalent material property of the honeycomb core

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

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

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

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

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

4 結論

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

(2)通過正交數值實驗設計(中心復合實驗設計)從蜂窩芯層的9個等效材料參數中篩選出3個重要參數作為有限元模型修正的輸入參量。

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