基于響應(yīng)面方法的印制電路板模型修正

徐飛，李傳日，姜同敏

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

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

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

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

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

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

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

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

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

4)測(cè)試結(jié)構(gòu)的振動(dòng)響應(yīng)特性，并利用有限元分析和試驗(yàn)結(jié)果構(gòu)造目標(biāo)函數(shù);

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

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

1.1 參數(shù)選擇

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

1.1.1 敏感性分析

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

1.1.2 重要度分析

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

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

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

1.2 試驗(yàn)設(shè)計(jì)方法

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

1.3 響應(yīng)面回歸分析

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

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

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

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

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

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

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

3)利用新群體更新設(shè)計(jì)點(diǎn)，計(jì)算目標(biāo)函數(shù)并檢驗(yàn)其收斂性，若滿(mǎn)足則修正完成，否則進(jìn)入下一步;

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

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

2 案例分析

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

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

2.1 模態(tài)試驗(yàn)

采用固定加速度傳感器，移動(dòng)力錘的方法對(duì)PCB進(jìn)行自由模態(tài)試驗(yàn)，共35個(gè)測(cè)試點(diǎn)(7×5)，試驗(yàn)示意圖和PCB測(cè)試點(diǎn)分布如圖3所示.

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

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

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

表2 試驗(yàn)振型的MAC值Table2 MAC values of test modes

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

2.2 有限元分析

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

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

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

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

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

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

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

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

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

圖7 響應(yīng)面示例Fig.7 Example of response surface

圖8 響應(yīng)面擬合度檢驗(yàn)Fig.8 Response surface fitness check

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

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

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

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

圖10 第4個(gè)目標(biāo)函數(shù)迭代過(guò)程Fig.10 Iteration process of fourth objective function

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

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

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

3 結(jié)論

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

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

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

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

本文雖然只研究了自由邊界條件下裸板的模型修正，但只要選擇好適當(dāng)?shù)男拚齾?shù)，本文提出的方法也可用于對(duì)真實(shí)使用條件下安裝了元器件的PCB進(jìn)行模型修正，這也是作者下一步將要研究的內(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］ 劉孝保，杜平安，夏漢良，等.一種面向動(dòng)態(tài)分析的PCB板等效建模方法［J］.儀器儀表學(xué)報(bào)，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).