基于響應面方法的印制電路板模型修正

徐飛，李傳日，姜同敏

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

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

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

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

1 基于響應面方法的模型修正技術

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

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

2)利用DOE構造采樣點并調用有限元模型計算輸出參數;

3)利用輸入參數和輸出參數構造響應面并對響應面進行回歸誤差分析;

4)測試結構的振動響應特性，并利用有限元分析和試驗結果構造目標函數;

5)在響應面模型的范圍內進行迭代分析，從而對目標函數進行優化，得到優化后的輸入參數和優化后的有限元模型.

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

1.1 參數選擇

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

1.1.1 敏感性分析

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

1.1.2 重要度分析

重要度分析用來確定各輸入參數對選定的輸出參數的影響，從而確定各輸入參數的重要程度.假設輸出變量與輸入變量之間存在線性或二次關系，利用R2準則來判定參數的重要度:

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

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

1.2 試驗設計方法

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

1.3 響應面回歸分析

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

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

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

1.4 目標函數構造與優化分析

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

1)創建初始群體;

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

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

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

圖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 模態試驗

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

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

圖3 試驗設置Fig.3 Test setup

圖4 速度頻響函數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個殼單元，邊界條件為自由狀態.有限元和試驗振型的MAC值如表4所示.表4中第1行表示試驗的前6階模態，第1列表示有限元計算的前6階模態.從表4可以看出，有限元和試驗的模態振型相關性較好(對角項均大于0.8)，且不同階振型的獨立性較好(非對角項均接近0).

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

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

式中，fai為有限元分析頻率;fei為試驗頻率，i為模態階數.

輸入參數的敏感性分析如圖5所示.

圖5 參數敏感性分析Fig.5 Parameter sensitivity analysis

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

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

圖6 參數重要度分析Fig.6 Parameter importance analysis

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

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

PCB的前3階模態頻率常常是最重要的，利用前3階共振頻率構造3個目標函數:

再利用前6階模態頻率殘差平方和構造第4個目標函數:

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

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

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

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

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

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

3 結論

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

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

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

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

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

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，Friswell 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，Friswell 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，Friswell 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，Fang 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］ 劉孝保，杜平安，夏漢良，等.一種面向動態分析的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).