徐飛，李傳日，姜同敏

(北京航空航天大學 可靠性與系統(tǒng)工程學院，北京100191)

有限元方法現(xiàn)被廣泛用于各領(lǐng)域的工程分析.然而，有限元方法是基于許多假設(shè)的，當這些假設(shè)與實際情況不符時會導致有限元分析與試驗結(jié)果之間存在較大誤差.這些誤差常常由以下幾方面的不確定性引起:結(jié)構(gòu)的材料屬性，幾何特征以及邊界條件等.模型修正技術(shù)是利用試驗結(jié)果改善有限元模型的質(zhì)量，盡可能減小有限元分析和試驗結(jié)果之間的誤差的技術(shù).模型修正技術(shù)的基本思想是通過分析輸出變量和輸入變量之間的關(guān)系，結(jié)合有限元分析和試驗結(jié)果構(gòu)造目標函數(shù)，利用優(yōu)化算法對有限元模型中不確定的輸入變量進行修正，從而修正有限元模型.早期的模型修正方法是直接修正法［1-5］.該方法的優(yōu)點是不需要進行迭代計算，計算量較小且不存在收斂性問題.這類方法能夠完全再現(xiàn)試驗結(jié)果，但修正后的質(zhì)量和剛度矩陣沒有確切的物理意義，且常常在工程上不可實現(xiàn).為了保留模型的物理意義，出現(xiàn)了迭代修正法［6-12］.迭代修正法需要在每次迭代過程中計算輸出變量對修正參數(shù)的敏感性，利用Taylor展開式計算參數(shù)的修正量.迭代修正法在每一次迭代過程中都要重新進行有限元計算，因此當模型較大時，會帶來耗時和收斂性差等問題.為了解決這樣的問題，近年來出現(xiàn)了基于替代模型的模型修正法.這種方法利用替代模型來替代原有限元模型，計算時不需要反復調(diào)用有限元模型進行迭代分析，因此計算速度更快，且收斂性較好.基于響應面的模型修正法是一種常用的替代模型法.該方法利用試驗設(shè)計法(DOE，Design of Experiment)對修正參數(shù)進行采樣并計算各采樣點的輸出，利用這些采樣點構(gòu)造響應面來近似輸出變量和輸入?yún)?shù)之間的關(guān)系，再利用有限元分析和試驗結(jié)果構(gòu)造目標函數(shù)，利用優(yōu)化算法對目標函數(shù)進行優(yōu)化從而得到參數(shù)修正量.基于響應面的模型修正方法常被用于大型結(jié)構(gòu)分析中，如橋梁，建筑等［13-15］.

在航空航天以及汽車等領(lǐng)域，有限元方法被廣泛用于印制電路板(PCB，Printed Circuit Board)的振動響應特性分析［16-21］.然而，PCB采用的是復合材料，其材料屬性在很大程度上取決于生產(chǎn)工藝.不同的PCB生產(chǎn)商生產(chǎn)出來的PCB材料屬性存在很大差異，且大部分情況下并不知道其確切值.雖然有部分文獻指出了PCB材料屬性的不確定性對其有限元分析結(jié)果有很大影響［22-24］，但目前很少有文獻給出詳細的PCB有限元模型的修正過程.

本文將給出利用基于響應面的模型修正技術(shù)對PCB有限元模型進行修正的過程，目標是利用自由邊界條件下的模態(tài)頻率對不確定的材料屬性(正交各向異性的彈性模量、主泊松比和剪切模量)進行估計.該修正過程首先將材料屬性作為輸入變量并對各輸入變量進行敏感性和重要度分析，確定出重要的輸入?yún)?shù);然后采用中心復合設(shè)計(CCD，Central Composite Design)進行試驗設(shè)計，對重要輸入?yún)?shù)進行采樣并計算;接著利用二次多項式響應面估計PCB的模態(tài)頻率和輸入?yún)?shù)的關(guān)系;最后利用有限元分析和模態(tài)試驗得到的模態(tài)頻率構(gòu)造多個目標函數(shù)，并采用遺傳算法進行優(yōu)化分析.本文給出一個案例對上述的修正過程進行了闡述.

1 基于響應面方法的模型修正技術(shù)

基于響應面方法的模型修正技術(shù)是一種基于響應面對結(jié)構(gòu)響應特性進行近似，進而對目標函數(shù)進行優(yōu)化分析，最終對不確定的輸入?yún)?shù)進行修正以提高有限元模型精度的技術(shù).基于響應面方法的模型修正技術(shù)流程圖如圖1所示.主要過程包括:

1)初步選擇優(yōu)化參數(shù)，確定其上下限，利用參數(shù)敏感性分析和重要度分析確定重要參數(shù)，從而確定最終的優(yōu)化參數(shù);

2)利用DOE構(gòu)造采樣點并調(diào)用有限元模型計算輸出參數(shù);

3)利用輸入?yún)?shù)和輸出參數(shù)構(gòu)造響應面并對響應面進行回歸誤差分析;

4)測試結(jié)構(gòu)的振動響應特性，并利用有限元分析和試驗結(jié)果構(gòu)造目標函數(shù);

5)在響應面模型的范圍內(nèi)進行迭代分析，從而對目標函數(shù)進行優(yōu)化，得到優(yōu)化后的輸入?yún)?shù)和優(yōu)化后的有限元模型.

圖1 基于響應面方法的模型修正技術(shù)流程圖Fig.1 Procedure of model updating technique based on response surface method

1.1 參數(shù)選擇

修正參數(shù)必須滿足兩個條件:①具有不確定性;②結(jié)構(gòu)響應對其敏感.當修正參數(shù)太多時，尤其是超過了可測量的結(jié)構(gòu)響應的數(shù)量時，優(yōu)化過程中將出現(xiàn)病態(tài)矩陣問題.因此，首先要對修正參數(shù)進行敏感性分析，確定少量的重要參數(shù)，提高計算精度和速度.

1.1.1 敏感性分析

敏感性分析能夠給出輸出變量對輸入變量的全局敏感性.正敏感性表示輸出變量隨著輸入變量的增大而增大，負敏感性表示輸出變量隨著輸入變量的增大而減小.本文采用Spearman秩排序法計算相關(guān)性系數(shù)以確定參數(shù)敏感性，該方法同時考慮了輸入?yún)?shù)和輸出參數(shù)的變化范圍.

1.1.2 重要度分析

重要度分析用來確定各輸入?yún)?shù)對選定的輸出參數(shù)的影響，從而確定各輸入?yún)?shù)的重要程度.假設(shè)輸出變量與輸入變量之間存在線性或二次關(guān)系，利用R2準則來判定參數(shù)的重要度:

式中，i為采樣點;N為總采樣點個數(shù);yi為在第i個采樣點的輸出變量值為在第i個采樣點的回歸值為所有yi的算術(shù)平均值.

R2越接近1，表明選定的輸入?yún)?shù)越重要，越接近0，表明該參數(shù)可從修正參數(shù)集中去掉，在有限元計算時作為常數(shù).

1.2 試驗設(shè)計方法

響應面是對采樣點的擬合，因此采樣點的選擇很大程度上決定了響應面的準確性和計算效率.采樣點少則導致響應面準確性差，采樣點多則導致計算效率低.在實際應用中，采樣點是利用DOE確定的.本文采用的DOE方法是常用的中心復合設(shè)計法，該方法能夠構(gòu)造出簡單的多項式型響應面.CCD采用正交表進行試驗設(shè)計以確定選定參數(shù)的采樣點.常用的CCD方法有3種:外切中心復合設(shè)計、嵌套中心復合設(shè)計以及面心立方設(shè)計.經(jīng)典的CCD方法包含3部分:①析因部分:一個立方體的2k頂點或這些析因點的一部分;②帶有參數(shù)α的2k個星點;③中心點.該方法是一種5水平部分析因設(shè)計方法:(－α，－1，0，1，+α)或(－1，－1/α，0，1/α，1).當星點處于設(shè)計空間每個面上時，CCD從5水平退化為3水平:(－1，0，1).

1.3 響應面回歸分析

響應面的類型對分析結(jié)果影響重大.不同的響應面可用于不同的應用情況.在與結(jié)構(gòu)動力學相關(guān)的模型修正中，多項式型響應面是常用的響應面類型，該類響應面計算簡單且精度較高.本文采用二次多項式響應面，計算公式如下:

其中，β0，βi，βij為回歸系數(shù);y 為響應面;x 為輸入變量.

采樣點必須大于等于多項式的項數(shù).當采樣點大于多項式的項數(shù)時，式(2)為超定方程，需要利用回歸技術(shù)對響應面進行擬合.在這種情況下，響應面通常不可能在每個樣本點得到精確解.最小二乘擬合常被用于構(gòu)造響應面.在利用響應面進行結(jié)果估計和模型修正之前，必須對響應面的擬合度進行檢驗.R2準則再次用來檢驗響應面擬合度.R2越接近于1，表明響應面擬合度越高，利用響應面估計得到的計算結(jié)果也就越準確.

1.4 目標函數(shù)構(gòu)造與優(yōu)化分析

在進行與結(jié)構(gòu)動力學相關(guān)的模型修正時，常常利用結(jié)構(gòu)的振動響應特性來構(gòu)造目標函數(shù).本文利用結(jié)構(gòu)的模態(tài)頻率來構(gòu)造目標函數(shù)，并利用多目標函數(shù)遺傳優(yōu)化算法(MOGA，Multiple Objective Genetic Algorithm)進行優(yōu)化分析.MOGA的流程圖如圖2所示，主要包括以下步驟:

1)創(chuàng)建初始群體;

2)利用交叉運算和變異運算創(chuàng)建新群體;

3)利用新群體更新設(shè)計點，計算目標函數(shù)并檢驗其收斂性，若滿足則修正完成，否則進入下一步;

4)檢查停止準則，若不滿足則回到第2)步繼續(xù)迭代，若滿足則停止，迭代失敗，重新定義目標函數(shù)，重新進行模型修正.

圖2 MOGA流程圖Fig.2 Procedure of MOGA

2 案例分析

案例分析對象為特制的菊花鏈電路PCB，尺寸為203 mm×140 mm×1.6 mm.初始的材料屬性如表1所示.

表1 PCB初始材料屬性Table1 Initial material properties of PCB

2.1 模態(tài)試驗

采用固定加速度傳感器，移動力錘的方法對PCB進行自由模態(tài)試驗，共35個測試點(7×5)，試驗示意圖和PCB測試點分布如圖3所示.

將測得的加速度頻響函數(shù)轉(zhuǎn)換為速度頻響函數(shù)，如圖4所示.利用模態(tài)指示函數(shù)(MIF，Mode Indicate Function)、穩(wěn)定圖(stabilization diagram)和全局時域復指數(shù)法提取共振頻率和振型.試驗振型的MAC(Modal Assurance Criterion)值如表2所示.從表2可以看出，測試點數(shù)量足夠多，振型之間相互獨立.試驗模態(tài)頻率的結(jié)果如表3所示.從表3可以看出，模型修正前的前6階頻率誤差最大為7.7%，最小為3.3%.

圖3 試驗設(shè)置Fig.3 Test setup

圖4 速度頻響函數(shù)Fig.4 Mobility

表2 試驗振型的MAC值Table2 MAC values of test modes

表3 修正前后前6階共振頻率對比Table3 Comparison of first six resonant frequencies before and after model updating

2.2 有限元分析

利用ANSYS對PCB進行有限元建模，采用200個殼單元，邊界條件為自由狀態(tài).有限元和試驗振型的MAC值如表4所示.表4中第1行表示試驗的前6階模態(tài)，第1列表示有限元計算的前6階模態(tài).從表4可以看出，有限元和試驗的模態(tài)振型相關(guān)性較好(對角項均大于0.8)，且不同階振型的獨立性較好(非對角項均接近0).

表4 有限元仿真和試驗振型的MAC值Table4 MAC values of finite element(FE)simulation and test modes

PCB的密度和幾何尺寸比較容易測量，且相對變化很小，因此本文將PCB的正交各向異性材料屬性(3個彈性模量，3個泊松比，3個剪切模量)作為初始輸入(修正)變量，前6階共振頻率和一個自定義變量作為輸出變量，自定義變量為前6階模態(tài)頻率的殘差平方和:

式中，fai為有限元分析頻率;fei為試驗頻率，i為模態(tài)階數(shù).

輸入?yún)?shù)的敏感性分析如圖5所示.

圖5 參數(shù)敏感性分析Fig.5 Parameter sensitivity analysis

由圖5可以看出，輸出變量只對Ex，Ey和Gxy敏感.利用這3個輸入變量進行重要度檢驗，結(jié)果如圖6所示.由圖6可以看出，當選擇這3個變量為重要參數(shù)時，前6階模態(tài)頻率和自定義輸出變量的R2均接近1，因此將這3個輸入?yún)?shù)作為最終的優(yōu)化參數(shù).利用DOE生成15個試驗點(采樣點)，生成算法為CCD.

利用2階多項式擬合出響應面，目標函數(shù)和Ex，Ey的響應面如圖7所示.利用R2準則進行響應面的擬合度檢驗，結(jié)果如圖8所示.由圖8可以看出，前6階共振頻率的R2均接近1，因此擬合度較高.

圖6 參數(shù)重要度分析Fig.6 Parameter importance analysis

圖7 響應面示例Fig.7 Example of response surface

圖8 響應面擬合度檢驗Fig.8 Response surface fitness check

PCB的前3階模態(tài)頻率常常是最重要的，利用前3階共振頻率構(gòu)造3個目標函數(shù):

再利用前6階模態(tài)頻率殘差平方和構(gòu)造第4個目標函數(shù):

利用MOGA對多目標函數(shù)迭代計算，目標是使得4個目標函數(shù)最小化.初始樣本數(shù)為100，每一次迭代生成100組新的樣本，迭代算法采用的是遺傳基因算法，每一次的變異概率(mutation probability)為 0.01，每一次交叉概率(crossover probability)為0.98.設(shè)置兩種收斂準則:最大允許Pareto比和收斂穩(wěn)定比.迭代過程如圖9～圖11所示.可以看出，經(jīng)過6次迭代后4個目標函數(shù)都達到了收斂值(第2個目標函數(shù)收斂性相對較差)，Pareto百分比達到穩(wěn)定值70%，Stability百分比達到穩(wěn)定值2%，收斂速度較快.

圖9 前3階共振頻率迭代過程Fig.9 Iteration processes of first three resonant frequencies

圖10 第4個目標函數(shù)迭代過程Fig.10 Iteration process of fourth objective function

圖11 兩種迭代判據(jù)下的收斂過程Fig.11 Convergence process with two different iterative criteria

修正前后PCB的材料屬性如表5所示.利用修正后的材料屬性得到修正后的模態(tài)頻率如表3所示.由表3可以看出，修正后PCB前6階模態(tài)頻率的誤差最大為3.2%，最小為0.

表5 修正前后材料屬性對比Table5 Comparison of material parameters before and after updating

3 結(jié)論

本文提出了利用DOE，二次多項式響應面，多目標函數(shù)和遺傳優(yōu)化算法對自由邊界條件下的PCB進行模型修正的方法，結(jié)論如下:

1)在自由邊界條件下，Ex，Ey和 Gxy對 PCB的模態(tài)頻率影響最大.

2)修正后PCB前6階模態(tài)頻率的誤差最大為3.2%，最小為0，達到了模型修正的目的.

3)多目標函數(shù)修正過程收斂速度較快，可利用已有商業(yè)有限元軟件直接進行分析，易于工程應用.

本文雖然只研究了自由邊界條件下裸板的模型修正，但只要選擇好適當?shù)男拚齾?shù)，本文提出的方法也可用于對真實使用條件下安裝了元器件的PCB進行模型修正，這也是作者下一步將要研究的內(nèi)容.

References)

［1］ Baruch M.Optimization procedure to correct stiffness and flexibility matrices using vibration tests［J］.AIAA Journal，1978，16(11):1208-1210.

［2］ Baruch M，Bar Itzhack I Y.Optimal weighted orthogonalization of measured modes［J］.AIAA Journal，1978，16(4):346-351.

［3］ Baruch M.Methods of reference basis for identification of linear dynamic structures［J］.AIAA Journal，1984，22(4):561-564.

［4］ Berman A.Optimalweighted orthogonalization ofmeasured modes-comment［J］.AIAA Journal，1979，17(8):927-928.

［5］ Berman A，Nagy E J.Improvement of a large analytical model using test data［J］.AIAA Journal，1983，21(8):1168-1173.

［6］ Fox R，Kapoor M，Rate of change of eigenvalues and eigenvectors［J］.AIAA Journal，1968，6(12)2426-2429.

［7］ Mottershead J E，F(xiàn)riswell M I.Model updating in structural dynamics:a survey［J］.Journal of Sound and Vibration，1993，167(2):347-375.

［8］ Friswell M I，Mottershead J E.Finite element model updating in structural dynamics［M］.Dordrecht:Kluwer Academic Publishers，1995:158-256.

［9］ Natke H G，Lallement G，Cottin N.Properties of various residuals within updating of mathematical models［J］.Inverse Problems in Engineering，1995，1(4)329-348.

［10］ Khodaparast H H，Mottershead J E，F(xiàn)riswell M I.Perturbation methods for the estimation of parameter variability in stochastic model updating［J］.Mechanical Systems and Signal Processing，2008，22(8):1751-1773.

［11］ Shahverdi H，Mares C，Wang W，et al.Clustering of parameter sensitivities:examples from a helicopter airframe model updating exercise［J］.Shock and Vibration，2009，16(1):75-88.

［12］ Mottershead J E，Link M，F(xiàn)riswell M I.The sensitivity method in finite element model updating:a tutorial［J］.Mechanical Systems and Signal Processing，2011，25(7):2275-2296.

［13］ Ren W X，F(xiàn)ang S E，Deng M Y.Response surface-based finiteelement-model updating using structural static responses［J］.Journal of Engineering Mechanics，2010，137(4):248-257.

［14］ Ren W X，Chen H B.Finite element model updating in structural dynamics by using the response surface method［J］.Engineering Structures，2010，32(8):2455-2465.

［15］ Zhou L R，Yan G R，Ou J P.Response surface method based on radial basis functions for modeling large-scale structures in model updating［J］.Computer-Aided Civil and Infrastructure Engineering，2013，28(3):210-226.

［16］ Pitarresi，J M.Modeling of printed circuit cards subject to vibration［C］//IEEE Proceedings of the Circuits and Systems Conference.Piscataway，NJ:IEEE，1990，3:2104-2107.

［17］ Pitarresi J M，Celetka D，Coldwel R，et al.The smeared properties approach to FE vibration modeling of printed circuit cards［J］.ASME Journal of Electronics Packaging，1991，113(3)250-257.

［18］ Pitarresi J M，Primavera A.Comparison of vibration modeling techniques for printed circuit cards［J］.Journal of Electronic Packaging，1992，114(4)378-383.

［19］ Lim G，Ong J，Penny J.Effect of edge and internal point support of a printed circuit board under vibration［J］.ASME Journal of Electronic Packaging，1999，121(2)122-126.

［20］ Deiter A，Baker M.Using analysis to design for drop or other shock environment［J］.Sound ＆ Vibration，2000，36(10):26-30.

［21］ Li R S.A methodology for fatigue prediction of electronic components under random vibration load［J］Transactions of the ASME，2001，123(4)394-400.

［22］ Amy R A，Aglietti G S，Richardson G.Sensitivity analysis of simplified printed circuit board finite element models［J］.Microelectronics reliability，2009，49(7):791-799.

［23］ Amy R A，Aglietti G S，Richardson G.Accuracy of simplified printed circuit board finite element models［J］.Microelectronics Reliability，2010，50(1):86-97.

［24］ 劉孝保，杜平安，夏漢良，等.一種面向動態(tài)分析的PCB板等效建模方法［J］.儀器儀表學報，2011，32(4):863-869.Liu X B，Du P A，Xia H L，et al.Dynamic property analysis-oriented PCB equivalent modeling method［J］.Chinese Journal of Scientific Instrument，2011，32(4):863-869(in Chinese).