孫衛(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)來進(jìn)行抑制。

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

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

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

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

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)過程，進(jìn)一步綜合考慮了軸向和剪切變形，推導(dǎo)出了精確的面內(nèi)泊松比[20]

(10)

對于剪切模量G23，可以通過對蜂窩芯單胞進(jìn)行有限元分析并對比不同高寬比的剪切模量，得到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]針對某飛機(jī)模型的前10階模態(tài)頻率建立了四階多項式響應(yīng)面模型。

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

1.2.1 實驗設(shè)計方法

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

對于確定性計算分析(包括有限元分析)不存在隨機(jī)誤差，而且隨機(jī)采樣可以最小化設(shè)計區(qū)域上的累計均方差[28]。目前最常用的隨機(jī)抽樣方法是拉丁超立方體設(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é)果，認(rèn)為在大多數(shù)情況下，RBF是在準(zhǔn)確度和魯棒性方面最可靠的方法，特別是對于高階非線性問題而言。

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)面模型進(jìn)行自適應(yīng)采樣的全局優(yōu)化算法。Jones等[30]最早提出了基于Kriging響應(yīng)面模型進(jìn)行自適應(yīng)LHD采樣進(jìn)行全局優(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)有采樣點間的距離容差，即當(dāng)新采樣點與現(xiàn)有采樣點的距離相對現(xiàn)有采樣點空間坐標(biāo)的百分比小于該容差限值時，優(yōu)化終止。

評估添加新樣本的可能最優(yōu)區(qū)域主要基于兩方面的考慮：(1)響應(yīng)面模型是否能夠預(yù)測到更好的設(shè)計；(2)在該區(qū)域內(nèi)的樣本數(shù)是否太少以至于響應(yīng)面模型在該區(qū)域的準(zhǔn)確性無法保證。因為這種優(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)特性。被測蜂窩板用軟橡膠帶懸掛以模擬自由邊界條件?？偣灿?6個敲擊點(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)的方式進(jìn)行蜂窩板建模，其中兩個面板使用殼單元，芯層使用實體單元。整個模型包含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)置信準(zhǔn)則(MAC)值均大于0.9，最大頻率誤差不超過3%，平均誤差為1.24%。說明有限元分析結(jié)果與實驗結(jié)果吻合較好。但對于第1階、第2階以及第5階模態(tài)，誤差在1%～3%之間，可期待通過模型修正進(jìn)一步提高其精度。

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

3 有限元模型修正

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

在進(jìn)行有限元模型修正之前，首先需要確定待優(yōu)化參數(shù)。對于蜂窩板蒙皮的材料參數(shù)、預(yù)埋件及膠層的質(zhì)量均可認(rèn)為是確定的?？紤]到蜂窩板加工工藝的誤差，應(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ìn)行計算，然后根據(jù)計算產(chǎn)生的樣本進(jì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和ρ)作為進(jìn)一步優(yōu)化的輸入?yún)?shù)，將有限元分析和實驗結(jié)果的前六階模態(tài)頻率的平均誤差作為最小化目標(biāo)函數(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)特征頻率，下標(biāo)i是模態(tài)階數(shù)。

使用徑向基函數(shù)響應(yīng)面模型進(jìn)行自適應(yīng)拉丁超立方采樣，進(jìn)行全局優(yōu)化。徑向基函數(shù)采用立方函數(shù)，每次迭代的采樣點數(shù)為10。優(yōu)化終止條件中新采樣點與現(xiàn)有采樣點的距離容差設(shè)為0.1%。圖7為整個迭代過程中的采樣點分布，其中不同顏色代表密度的大小(范圍為20～35kg/m3)，圓的直徑尺寸代表頻率誤差的大小(范圍為0.88%～22.24%)，從圖7中可以看到當(dāng)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為當(dāng)密度為23.58 kg/m3時，G23和G31相對模態(tài)頻率平均誤差的響應(yīng)面模型。圖9中可以看出，當(dāng)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ù)對蜂窩芯的兩個面外剪切模量及密度進(jìn)行了參數(shù)優(yōu)化，修正后的有限元模型的前六階模態(tài)的平均頻率誤差達(dá)到0.88%，最大誤差不超過2.5%，精度可以滿足航天器微振動分析的要求。