局部涂敷硬涂層薄板有限元建模及涂敷位置優化

劉 蓉， 孫 偉

(東北大學 機械工程與自動化學院，沈陽 110819)

硬涂層是金屬基、陶瓷基或者兩者混合制成的涂層材料，主要用于熱障、耐摩擦和抗腐蝕涂層。近年來的研究[1-3]發現硬涂層還具有阻尼減振效果，表現為：在結構表面噴涂硬涂層可顯著降低基體結構的共振應力幅度。進一步的研究發現，硬涂層之所以能減振是源于硬涂層顆粒之間的內摩擦，Tassini[4]，Torvik[5]，Al-Rub[6]等分別創建了微觀材料學表征模型來解釋硬涂層的減振機理。在實際應用中，在結構表面全部涂敷硬涂層是不現實的。一方面，大量結構不允許在減振過程中顯著增加或改變其自身的質量、體積和固有頻率，另一方面，由于功能限制，通常僅可以在結構的一部分表面上沉積硬涂層。可見，為了有效實施硬涂層減振，除了需要系統地研究涂敷硬涂層復合結構動力學建模(尤其是局部涂敷硬涂層復合結構的建模)，還要研究硬涂層復合結構的阻尼優化問題。

關于硬涂層復合結構動力學建模與分析，學者們已開展了大量研究。Yang等[7]推導了硬涂層復合板的控制方程，并且基于Lindstedt-Poincaré攝動法求解其固有特性。Li等[8]針對硬涂層懸臂薄板結構，研究其非線性振動機理并且采用有限元迭代法計算其固有頻率和振動響應。Chen等[9]研發了一種解析方法以計算葉片上涂敷硬涂層的整體葉盤的自由振動特性及阻尼效應。上述關于硬涂層復合結構的建模與分析主要針對的全涂敷結構，即在可涂層區域完全涂敷硬涂層。硬涂層復合結構的阻尼優化需要局部涂敷模型，而有關局部涂敷涂層復合結構的建模與分析研究還不充分。

作為一項新興的減振技術，有關硬涂層復合結構阻尼優化的研究并不多見，但是與之相類似的黏彈性阻尼結構優化已開展了大量的研究。Lumsdaine等[10-12]針對梁和板型結構以峰值位移最小化(或者系統損耗因子最大化)為目標進行形狀優化。Chen等[13]以系統的結構阻尼作為主要性能指標，對約束阻尼材料的布局進行了拓撲優化。Hou等[14]采用遺傳算法確定約束阻尼梁結構約束層和黏彈性層的厚度以及黏彈性層的剪切模量的最優值。Zheng等[15]以圓柱殼振動響應最小化為目標，采用基于罰函數法的基因優化對被動約束阻尼層進行布局優化。

上述關于黏彈性阻尼層貼敷位置的優化獲得的優化結果通常是離散的。硬涂層的制備工藝主要是等離子噴涂和物理氣相沉積[16]，而用上述兩種工藝很難按離散的優化方案在基體上沉積硬涂層。為此，結合涂層工藝，這里對硬涂層復合結構的阻尼優化定義為：在固定硬涂層形狀的基礎上，尋找最佳的涂層位置。這樣，只需在基體上一個指定的連續區域(涂層的形狀及面積是固定的)沉積涂層，從而使硬涂層阻尼優化在工程上更加便于實施。本文面向局部涂敷硬涂層薄板結構，研究以遺傳算法求解該優化問題的方法。推導了局部涂敷硬涂層薄板結構的有限元方程，確定了以優化涂層位置為目標的優化模型，最后以涂敷NiCrAlCoY+YSZ硬涂層材料的懸臂板為例進行了實例研究。

1 局部涂敷硬涂層薄板有限元建模及阻尼預估

1.1 局部涂敷硬涂層薄板結構有限元建模

如圖1(a)所示的處于懸臂狀態的局部涂敷硬涂層的薄板結構，圖1(b)為薄板涂層區域部分橫截面，Hc和Hs分別為硬涂層及基體的厚度，d為復合結構中性層到基體中性層的距離，Hsc為基體和硬涂層中面的距離。這里采用等效單層理論對其進行有限元建模。由于是局部涂敷硬涂層，將該硬涂層板劃分為若干個單元后，整個系統會出現兩種單元：涂層復合結構單元(同時包含基體和涂層)和基體單元。為了方便，這里用4節點板單元同時模擬上述兩種單元，只是針對不同的單元類型賦予不同的材料參數。單元的具體結構見圖2，每節點具有三個自由度，具體為一個位移自由度w(撓度)、兩個轉動自由度θx和θy， 單元的長度及寬度分別為a和b。

圖1 局部涂敷硬涂層的薄板結構Fig.1 A partially-covered hard-coating cantilever thin plate

圖2 四邊形板單元Fig.2 Quadrilateral plate element

設δe單元節點位移向量，可表示為

δe={w1θx1θy1…w4θx4θy4}T

(1)

按有限元理論可確定單元的形函數矩陣N，進一步，按照幾何方程可推導出單元的應變矩陣B，表示為

B=[B1,B2,B3,B4]

(2)

其中

(3)

式中：Bi是個3×3矩陣， 因而應變矩陣B是一個3×12的矩陣。

硬涂層的內阻尼雖然大于金屬但遠小于黏彈性阻尼材料(硬涂層的損耗因子一般是金屬的2倍～10倍)，為了精確描述涂層結構的振動特性，金屬基體的內阻尼也應考慮。這里硬涂層及金屬基體的材料參數均用復模量表示，具體為

(4a)

(4b)

(5)

式中：vs為基體的泊松比。

對于包含硬涂層的復合單元，彈性矩陣為基體和硬涂層的耦合，因而復合層結構的等效彈性矩陣可表示為

(6)

(7)

硬涂層對復合彈性矩陣的貢獻為

(8)

式中：vc為硬涂層的泊松比，而復合結構中性層到基體中性層的距離d可表達為

d=EcHcHsc/(EsHs+EcHc)

(9)

在確定單元的應變矩陣與彈性矩陣后，可按虛位移原理得到單元剛度矩陣。對于不含硬涂層的單元，剛度矩陣可表示為

(10)

對于包含硬涂層的單元，剛度矩陣可表示為

(11)

與剛度矩陣相對應，局部涂層薄板結構單元質量矩陣的求解式也有兩種。對于只有基體的單元，質量矩陣的求解式為

(12)

對于同時包含基體及硬涂層的復合單元，質量矩陣的求解式為

(13)

式中：ρs,ρc分別為基體及硬涂層的密度。

引入邊界條件，并對單元剛度與質量矩陣進行組集，從而求得總剛度矩陣K*和質量矩陣M，最終可獲得如下頻域運動方程

(K*-ω2M)x=F

(14)

式中：ω為激振頻率；x為結構總的位移響應向量；F為激振力向量。

1.2 基于修正模態應變能法的阻尼預估

準確求解復合結構的模態損耗因子是實施阻尼優化的基礎。眾所周知，經典的模態應變能法[17]是預估復合結構模態損耗因子的經典方法。參照式(14)，依據經典模態應變能法用于求解復合結構模態振型的特征方程可表示為

(15)

式中：KR為復剛度矩陣K*的實部，ωr,φr分別是r階固有頻率和模態陣型。最終，按照經典模態應變能法硬涂層復合結構的阻尼預估公式可表示為

(16)

式中：KIc為硬涂層對復剛度矩陣K*虛部(KI)的貢獻。

上述經典模態應變能法通常適用于金屬基體表面貼敷黏彈性阻尼材料的復合結構。因為黏彈性材料的損耗因子遠大于金屬(大約是金屬的100倍～1 000倍以上)，因而在計算中不用考慮金屬的損耗因子。但是，該經典方法并不能直接用于硬涂層復合結構的阻尼預估，原因有二，一是式(16)中只考慮了硬涂層對復剛度矩陣虛部的貢獻，由于硬涂層的損耗因子與金屬較為接近，因而對于該復合結構金屬基體的阻尼也不應忽略；二是求解復合結構模態振型時僅考慮了復剛度矩陣的實部。為此，需要對經典的模態應變能法進行修正，以確定適用于硬涂層復合結構的阻尼預估。

這里，將求解復合結構模態振型的特征方程式(15)修正為

(17)

(18)

式中： trace()是求矩陣的跡的函數, 即矩陣的對角線元素之和。

最終，基于修正模態應變能法，硬涂層復合結構阻尼性能預估公式變為

(19)

式(19)中同時考慮基體及硬涂層損耗因子的影響，而且考慮復剛度矩陣的虛部對模態振型的貢獻，因而其阻尼性能的預估精度應高于經典的模態應變能法。由于硬涂層自身的損耗因子大約在10-2或10-3這一數量級，且相對于金屬基體厚度很薄，因而硬涂層復合結構仍滿足小阻尼特性，但獲得的各階模態損耗因子仍大于金屬的損耗因子。

2 阻尼優化

2.1 優化設計模型

以下以薄板的有限元模型為例描述為實現上述優化，所定義的優化模型。薄板的有限元模型見圖3，這里共分了256個單元，黑色區域表示涂層單元。

圖3 局部涂敷硬涂層的薄板結構有限元模型Fig.3 The finite element model of the partially-covered hard-coating cantilever thin plate

這里將設計變量定義為矩形區左下角單元的位置p，可以用單元編號來描述。由于涂層形狀是固定的，因而對于一個分網完畢的結構，位置p并不能等于任意單元編號，其可以等于的編號可組成一個向量P， 表達為

p∈P=[p1,p2,…,pi,…,pN]

(20)

式(20)中N為位置p可達的總數，其可以由所規劃的涂層形狀來確定。進一步，一旦指定了位置p，可以由預設的涂層形狀快速確定整個涂層區域，例如，圖3中p=33，整個涂層區規劃為9×9，則可按照單元編號規律確定所有其他的涂層單元編號。

優化的目標函數可描述為：指定連續的涂層區后系統的某一階或者若干階具有較大的模態損耗因子。因而最終的優化模型可定義為

(21)

式中：ar為對應優化目標的加權系數。

2.2 優化算法

這里采用多種群遺傳算法(Multiple Population GA,MPGA)[18]來求解式(21)所描述的優化模型。該算法中多個種群使用同一目標函數，種群之間定期進行信息交換，各種群的交叉率和變異率取不同的固定值以搜索不同解空間的最優解。

MPGA方法優化步驟具體敘述如下：

(1) 編碼。采用二進制編碼方式，將位置P可達的總數N編碼成S位二進制串。例如，N=123，則可以編碼成22位二進制串,如圖4所示。

圖4 一個染色體二進制串Fig.4 A chromosome

(2) 種群初始化。隨機產生n個初始種群，每個種群含有m個染色體，每個染色體對應一個方案解。

(3) 適應度評價。根據目標函數f(p)合理定義適應度函數，并計算各初始種群的個體適應度值。

(4) 根據適應度值的大小選擇下一代的種群，適應度值表明個體的優劣，適應度越大，被選中的概率就越高，最合適的個體生存下來并進行復制，不合適的個體則被淘汰。

(5) 對各種群執行相應的交叉和變異操作，從而產生下一代種群。

(6) 各種群相互獨立，協同進化，相互之間通過移民算子聯系。在進化的每一代，通過人工選擇算子選出各種群中最優個體放入精華種群加以保存，最終采用最優個體最少保持代數作為終止判據，結果輸出最優個體，即最佳涂層左下角單元編號。

3 研究實例

這里以單面局部涂敷NiCoCrAlY+YSZ硬涂層材料的懸臂薄板為例實施硬涂層涂敷位置的優化。

3.1 研究對象

所研究的鈦板試件見圖5，相關的幾何及材料參數見表1。如“2.1”節所述，薄板共劃分了16×16個單元，整個涂層區規劃為9×9。需要說明的是表1中所列的尺寸為薄板參與振動分析的有效尺寸，而懸臂板夾持區的長度為30 mm。

圖5 鈦板實驗件(mm)Fig.5 Titanium plate specimen (mm)

表1 復合板相關尺寸及材料參數

3.2 局部涂敷硬涂層薄板振動測試與分析

利用錘擊法對涂層懸臂鈦板進行測試，圖6為局部涂敷硬涂層的懸臂板固有頻率測試的現場圖。所用到的儀器設備包括：PCB模態力錘，Polytec激光多普勒測振儀，LMS 16通道便攜式數據采集前端控制器和LMS.Testlab筆記本工作站等。由測試獲得的頻響函數可獲得涂層鈦板的固有頻率及阻尼比。利用Matlab按照第2部分所描述的原理自編有限元程序計算該局部涂層板的固有頻率，并將獲得的結果與實驗和基于ANSYS的計算結果進行比較。圖7為基于ANSYS所創建的有限元模型(在單元數量上與自編有限元程序同一規模)，表2為具體的實驗和分析計算結果。

圖6 局部涂敷硬涂層的懸臂板測試的現場圖Fig.6 Actual photograph of hammering test

圖7 局部涂敷硬涂層懸臂薄板有限元模型Fig.7 The finite element model of the partially-covered hard-coating thin plate

對比用Matlab自編有限元和基于ANSYS的計算結果可以看出，固有頻率計算的最大偏差為0.28 %；而基于Matlab計算獲得的固有頻率與實驗結果的最大偏差為1.79%，因而，可以證明本文自編有限元算法的合理性。需要說明的是，這里并沒有比較自編有限元計算和實驗獲得的硬涂層薄板阻尼參數，因為，圖6中實驗獲得的阻尼不僅包含了硬涂層及基體的材料阻尼，還包含夾持區及空氣中的阻尼，兩者在數值應該不一致。

表2 MATLAB自編有限元，基于ANSYS軟件和實驗獲得的

3.3 硬涂層涂敷位置優化

基于第2節建立的優化模型及算法對局部涂敷硬涂層薄板進行阻尼優化。優化是以獲得局部涂敷硬涂層板最大模態損耗因子為目標，包括單階次及多階次優化。在所有優化計算中，采用如下設置：種群數目為3，每個種群個體數為10，在[0.7 0.9]范圍內隨機產生交叉概率，在[0.001 0.05]范圍內隨機產生變異概率。另外，設定當目標函數值滿足10代不發生變化時遺傳算法終止。

3.3.1 單階次優化

這里，分別以復合結構前5階模態損耗因子最大為優化目標實施單階次優化。具體地，在針對r階優化中，將式(21)優化模型中的加權系數ar設置為1，其余加權系數置0，進而實施優化計算。在Intel Core i7-6700 CPU的計算機平臺上單階次優化及后續的多階次優化所耗費的時間大約25～30 min。相應的優化結果見圖8(圖中數字是指涂層左下角單元編號)，在表3中給出了對應各階次優化的涂層方案復合結構的前5階模態損耗因子。從表3中可以發現，在上述5種針對單階次的阻尼優化中，目標階次的模態損耗因子一般都較大。

表3 針對單階次優化后獲得的相應各階模態損耗因子(10-3)

圖8 針對單階次優化時涂層位置分布Fig.8 The location of the coating when performing optimization for single order

3.3.2 多階次優化

這里以同時優化1、3、5階和3、4、5階兩種情況為例描述面向硬涂層薄板多階次優化的實現過程。在針對1、3、5階進行優化時，相應地取各階次的加權系數分別為0.7、0.1、0.2；而對3、4、5階進行優化時，加權系數同樣也取為0.7、0.1、0.2。某階次的加權系數大，意味著這一階的模態損耗因子最大化在優化過程中被重點考慮。最終的優化結果見圖9，而對應優化結果的硬涂層板各階次模態損耗因子見表4。從表4可以看出，多階次優化后復合板所有目標階次對應的模態損耗因子較大，這也進一步說明了所研發的優化模型及算法的合理性。

圖9 針對多階次優化時涂層位置分布Fig.9 The location of the coating when performing optimization for multi orders

表4 針對多階次優化后獲得的相應各階模態損耗因子(10-3)

3.4 實驗驗證

為說明上述優化結果的合理性，這里以面向第1階的優化為例，采用實驗方法對優化結果進行驗證。圖10為具有不同涂層方案的3塊板，其中圖10(a)為對應第1階的優化方案，圖10(b)和(c)為在非第1階優化區域隨機選取位置涂敷硬涂層的阻尼方案。采用錘擊法激勵并用半功率帶寬法辨識復合結構的阻尼，相應的結果見表5。從表5可以看出，對應優化方案獲得的第1階模態阻尼比明顯大于其他兩種非優化方案的，因而可以證明針對第1階的優化結果是合理的。客觀上，也能說明本文研發的優化模型及算法也是合理的。

圖10 具有不同阻尼方案的3塊硬涂層板Fig.10 Three hard-coating plates with different damping schemes

表5 測試獲得的3塊涂層后鈦板的模態阻尼比

4 結 論

面向局部涂敷硬涂層的薄板結構，本文創建了該復合結構的有限元模型，并研究了局部涂層薄板阻尼優化的方法，得出以下結論：

(1) 采用等效單層理論對局部涂敷硬涂層復合結構進行有限元建模，采用四邊形板單元推導出了復合板的剛度矩陣和質量矩陣。

(2) 這里用修正模態應變能法對硬涂層復合結構進行阻尼預估，在該方法中同時考慮了基體及硬涂層損耗因子的影響，而且考慮復剛度矩陣的虛部對模態振型的貢獻，因而其阻尼性能的預估精度應高于經典的模態應變能法。

(3) 將涂層塊的左下角單元編號定義為設計變量，以獲得單階次或者多階次最大模態損耗因子為目標確定本文的硬涂層薄板阻尼優化模型并提出利用多種群遺傳算法求解該優化問題。實踐顯示，利用上述優化模型，無論是執行單階次還是多階次優化，均可以發現目標階次阻尼最優對應的涂層位置。