基于模型修正的螺栓連接板結構動力學建模方法

郭 寧，黃 慶,2，曹善成，徐 超

(1. 西北工業大學航天學院，西安 710072；2. 湖北航天飛行器研究所，武漢 430040)

一般情況下，影響飛行器結構動力學建模準確性的一個主要因素是其連接部位的仿真預示[1]。螺栓連接板結構是飛行器結構中最常用最重要的一種可拆卸式連接結構形式。因此，準確建立螺栓連接板結構的動力學模型，是準確獲取飛行器結構動力學特性和響應的前提之一。

目前，針對螺栓連接部件動力學建模的主要方法有：線性建模方法、非線性建模方法及薄層單元方法。線性方法是將實際的非線性連接結構近似成線性對象進行動力學分析，主要有三種方法：全固連建模方法[2]、節點約束方法[3]和用桿單元、彈簧單元和梁單元等模擬結構連接形式[4?5]。線性建模方法具有工程操作簡單等優點，但是不能準確地表征結構連接部位的剛度和阻尼特性，導致動響應分析結果偏差較大。非線性方法是通過假設或采用試驗途徑來建立連接界面的唯象力學模型。目前應用較多的界面力學模型有：Iwan模型[6?7]、LuGre毛刷模型[8]和Valanis模型[9]等。上述非線性方法存在非線性機理描述有限、非線性模型參數獲取難和計算耗費大等問題，還未有效地應用于工程結構分析中。薄層單元法最初是用于建立巖石與巖石界面之間的接觸模型。目前已被廣泛應用于連接結構建模中，Mayer和Gaul[10]的研究表明薄層單元能夠有效地模擬螺栓連接接觸面的力學特性。翟學等[11]研究了基于薄層單元法的螺栓連接結構建模和分析方法，并應用于渦扇發動機的建模。Iranzad和Ahmadian[12]使用實驗數據(固有頻率、頻響函數)來確定薄層單元的材料本構參數。姚星宇等[13]給出了薄層單元法在航空發動機機匣螺栓連接結構的建模方法，并研究了薄層單元材料參數對螺栓連接剛度的影響規律。王攀和臧朝平[14]基于螺栓連接超模型的剛度理論、赫茲接觸理論以及M-B分形模型，考慮了螺栓的數量和法蘭邊的接觸因素，建立了基于分區域薄層單元的螺栓連接簡化建模方法。姜東等[15? 16]建立了基于各向同性本構關系的薄層單元理論，提出了螺栓連接結構接觸面不確定性參數的識別方法。從上述的研究可以看出，采用薄層單元法建立連接結構動力學模型時薄層單元的本構參數的確定至關重要。

本文結合薄層單元理論，針對螺栓連接板結構提出了基于模型修正的動力學建模方法，其中采用薄層單元的力學行為等效模擬螺栓連接接觸界面的力學行為，建立用于模型修正的初始動力學有限元模型，將不確定的薄層單元材料本構參數作為修正變量，通過模態參數構建目標函數，提出采用水循環算法(Water Cycle Algorithm, WCA)來修正薄層單元的不確定本構材料參數。

1 基于薄層單元法的螺栓連接板結構動力學建模

圖1所示為工程中典型的螺栓連接板結構，首先采用板(殼)單元對被連接板部件進行動力學建模；繼而采用一定厚度的8節點六面體實體單元(即薄層單元)來等效模擬螺栓連接力學行為；最后在薄層單元與對應接觸的殼單元之間建立綁定約束關系，建立綁定約束時，需要設置綁定容差，容差應略大于兩個被連接板中面之間的距離，以保證薄層單元與板(殼)單元之間能夠有效的約束；螺栓的質量采用集中質量單元進行模擬。具體如圖1所示。

圖 1 螺栓連接板結構動力學建模Fig.1 Dynamic modeling of bolted joint plate

根據上述建模過程可知，螺栓連接板結構的連接界面剛度是通過薄層單元的剛度來等效的，建模過程中需要確定薄層單元的厚度以及其材料本構。薄層單元的厚度對分析結果影響較大，厚度過大則單元會有6個應變分量，難以準確體現接觸界面的力學特征；厚度過小則會導致單元雅可比矩陣行列式值趨向于零，造成矩陣病態，求逆困難，無法準確計算其位移-應變關系[15]。本文采用文獻[17]中給出的薄層單元厚度計算方法，即：

式中：l1和l2分別為薄層單元的長度和寬度；t為薄層單元厚度；R的取值范圍為10～100。

由于薄層單元的厚度遠小于其長度和寬度，因此，可忽略薄層單元的面內應變分量( εx,εy,γxy)。并且假設連接界面的法向和切向接觸性能相互獨立，且2個切向的接觸性能相同，則表征界面接觸性能的薄層單元本構方程可表示為[15]：

式中，En和Gt分別是薄層單元的法向彈性常數和切向剪切常數。若接觸面的法向與切向接觸性能相互耦合，可在本構關系式(2)中通過引入耦合項來實現。

類比方程式(2)，在有限元建模中，可采用各向同性材料的本構關系近似替代薄層單元本構關系，具體為[15]：

式中：G為材料剪切模量；λ為拉梅常數。

式中：E為材料彈性模量；μ為材料泊松比。

在實際情況中，由于螺栓預應力的作用，螺栓連接界面的剛度分布并非一致，螺栓附近區域的界面連接剛度高于遠離螺栓區域，因此，將薄層單元分為兩部分，采用不同的材料本構參數進行建模。分區大小由螺栓預應力的作用區域大小決定，如圖2所示，圖中 α角的取值范圍為25°～33°[18]，本文選擇 α=30?。為了進一步簡化有限元建模，將圓形的螺栓預應力作用區域等效為正方形區域來進行有限元網格離散，但需要保證作用區域面積相同。

圖 2 螺栓連接界面應力分布示意圖Fig.2 Stress distribution at bolted joint interface

2 基于水循環算法的螺栓連接板結構動力學模型修正

根據第1節的分區薄層單元法建立的螺栓連接板結構動力學模型可知，需要確定的建模參數包括：螺栓預應力作用區域的薄層單元材料參數E1、μ1以及非作用區域的薄層單元材料參數E2、μ2。

2.1 目標函數的構建

本文利用固有頻率的誤差和模態置信準則殘差的加權和來構建模型修正的目標函數：

式中：N表示實驗測試模態階數；p={p1,p2,···,pm}T為待修正變量向量，其中，下標m為修正變量個數；ωiE和ωiA分別為實驗和仿真的第i階固有頻率；α為加權系數； MACi為第i階模態的模態置信度。

2.2 水循環優化算法

從式(6)可以看出，當目標函數 Π(p)取得最小值時，即可得到模型修正問題的解，即薄層單元的待確定參數。因此需要將式(6)轉化為如下的最小值問題：

本文選擇水循環優化算法來執行式(9)所示的最小值優化問題。水循環優化算法是Eskandar等[19]于2012年提出的一種新穎的優化算法，其理論受啟發于大自然，是通過觀察大自然水循環過程中水的流動和轉換而提出的。目前，水循環算法已經被廣泛地應用于各類實際工程優化問題中，如Cho和Kim[20]將水循環算法應用于混合電磁鐵的優化設計方法中，通過實驗驗證了該算法的性能優于常規優化算法；Ghaffarzadeh[21]提出采用水循環算法來開發電力系統穩定器，結果表明，采用水循環算法可以在較短時間內檢測出最優參數，進一步提高了電力系統的穩定性。然而，據本文作者所知，水循環算法尚未被應用到解決結構動力學模型修正的問題中。

執行水循環算法，首先要隨機生成一個初始降雨層(類似于遺傳算法的種群或粒子群算法中的粒子群)，其包含的個體稱為雨滴(Raindrop)，代表優化問題中的一組解(類似于遺傳算法的染色體或者粒子群算法中的粒子)。初始降雨層中的每個雨滴Ri是在給定范圍內隨機產生的向量，對于本文中的m維修正變量問題，其表達式為：

式中： δi∈U(0,1) 為服從均勻分布的隨機數；n為降雨層中的雨滴個數。

其次，將所有雨滴根據其目標函數值(即式(6))的大小劃分為三個等級，其中最好的雨滴(即當前目標函數值最小的雨滴)形成海洋，較好的一些形成河流，其余的雨滴形成流向河流和海洋的溪流：

式中，下標NSr為河流數和海洋數之和，且海洋數為1。

假定所有溪流都會匯入河流或者海洋，由于河流和海洋的流量不同，流入每條河流和海洋的溪流數量NS就不同，具體為：

式中，NStream為溪流的數目。

然后，水循環算法進入匯流階段。其中：一部分溪流先匯入河流，再匯入海洋；另一部分溪流則直接匯入海洋。具體匯流的表達式如下：

為了防止算法陷入局部最優，通過蒸發和降雨來擴大算法的搜索空間。降雨過程可根據降雨位置的不同分為兩種：一種是在溪流或河流附近降雨，即隨機產生溪流或河流的新個體，增大算法的搜索域；另一種是在海洋附近降雨，即在海洋附近隨機產生新個體，使算法在當前最優解附近尋找更優解。可以根據以下兩個條件來判斷是否滿足降雨條件：

式中，dmax為人為設置的閾值，通常在算法初始階段取dmax=10?16，且隨著迭代地進行自適應減小，具體為：

式中，M為最大迭代次數。

如果滿足式(16)，則在整個搜索空間內進行降雨，具體為：

如果滿足公式(17)，則在海洋附近區域進行降雨，具體為：

式中：λ為服從均勻分布的隨機數向量；ν為在海洋附近的搜索范圍的大小，通常取ν=0.1。

同其他啟發式智能優化算法一樣，本文選用最大迭代次數作為水循環優化算法的終止條件。

綜上，結合薄層單元法和模型修正方法的螺栓連接板結構動力學建模流程如圖3所示，且具體實施步驟如下：

第一步，建立部件結構的動力學模型，并對其建模參數進行修正；

第二步，采用薄層單元法建立螺栓連接板結構的初始有限元模型，并進行模態分析，提取仿真模態參數；進行螺栓連接板結構的實驗模態分析，提取相應的實驗模態參數；

第三步，利用實驗模態參數與對應仿真模態參數，構造模型修正目標函數；

第四步，選擇修正變量，給定修正范圍；

第五步，設置水循環算法的初始參數，如初始雨滴數目n，海洋加河流的數目NSr，溪流數目參數NStream，蒸發條件判斷數dmax，最大迭代次數M等；

第六步，產生初始種群，根據目標函數值劃分海洋、河流和溪流，計算匯流數目，進行相應的匯流；

第七步，判斷是否滿足蒸發條件，如果滿足，則執行降雨過程；

第八步，判斷是否滿足算法終止條件，如果滿足，則結束算法并輸出最優解；如果不滿足，則返回到匯流過程，繼續執行算法，直到滿足終止條件為止；

第九步，輸出算法最優解，即為建模所需的薄層單元的材料參數。

圖 3 基于模型修正的螺栓連接板結構動力學建模流程圖Fig.3 Flow chart of dynamic modeling of bolted plate structure based on model updating

3 實驗驗證

為了驗證本文所提的螺栓連接板結構動力學建模方法的可行性和實用性，設計了如圖4所示的兩端固支五螺栓搭接板模型，其中，單板的尺寸為240 mm×120 mm×3 mm(去除邊界加持區域)；搭接區域長度為40 mm。搭接板材料為鋁合金，采用5個M6不銹鋼螺栓進行連接，且擰緊力矩為5.2 N?m。

圖 4 五螺栓搭接板結構Fig.4 Five-bolt lap plate

實驗結構與測試設備如圖5所示，包括信號采集及模態識別系統(LMS SCADAS Mobile)、PCB力錘(靈敏度為2.267 mV/N)、PCB加速度傳感器(靈敏度為100.5 mV/g)和STANLEY扭矩扳手。實驗中采用“移動激勵點法”進行模態測試。為了抑制由于信號截斷引起的泄漏，實驗測試中對輸入信號(即激勵信號)添加力窗，并對輸出信號(即響應信號)添加指數窗。為了一定程度上消除測試噪聲的影響，實驗中在每個測點進行3次測量，并對3次測試響應數據進行平均。本實驗所關心的結構的頻帶范圍為1 Hz～620 Hz，包含了結構前6階模態，實驗中設定振動信號的采樣頻率為1280 Hz(大于最大測試頻率的2倍)，頻率分辨率為0.31 Hz。實驗獲得的結構前六階模態結果見表1。

圖 5 模態試驗現場布置圖Fig.5 Layout of modal test

首先對部件結構(即板1和板2)的模型參數進行修正。通過對實驗件進行測量，得到板1和板2的尺寸分別為240 mm×120 mm×2.96 mm和240 mm×120 mm×2.95 mm，重量分別為284 g和282 g，可計算得到結構的材料密度為2717.7×10?12t/mm3。所以僅需要對其材料彈性模量E和泊松比μ進行修正。選取板1為修正結構，進行模態實驗，邊界條件為自由-自由，且采用4邊形殼單元對單板結構進行有限元建模，修正后的材料參數見表2。

表 1 五螺栓搭接板前六階實驗模態Table 1 First six experimental modes of five-bolt lap plate

根據第2節中所述的螺栓連接結構的動力學建模方法，建立了如圖6(a)所示的有限元模型，薄層單元采用六面體8節點實體單元，其中白色單元區域為螺栓預應力作用區域，紅色單元區域為非螺栓預應力作用區域，如圖6(b)所示。對結構進行有限元離散時選擇的最大單元邊長為4.5 mm，根據式(1)，我們選取得薄層單元厚度為0.1 mm。

圖 6 五螺栓搭接板有限元模型Fig.6 Finite element model of five-bolt lap plate

定義設計變量集p={E1,μ1,E2,μ2}，每個設計變量相對于固有頻率的靈敏度如圖7所示。

從圖7可以看出，泊松比μ1和μ2對低階模態不敏感，因此，選擇敏感性較強的設計變量E1和E2作為模型修正的修正變量，并分別給定修正區間為E1∈[0.00110] 和E2∈[0.00110]；建模中設定μ1=μ2=0.3。

采用第2節建立的模型修正方法對修正變量進行修正，其中水循環算法的具體參數設定見表3。

圖 7 設計變量對于固有頻率靈敏度分析Fig.7 Sensitivity of natural frequencies versus design parameters

表 3 水循環算法的參數設定Table 3 Water cycle algorithm parameters

將修正后的參數輸入有限元模型(即E1=0.001 MPa，E2=0.2957 MPa)，進行模態分析，并將參數修正后的模態分析結果與實驗結果進行對比，詳見表4。

表 4 修正后模態分析結果Table 4 Modal analysis results after updating

從表4可以看出，采用本文提出的動力學建模方法建立的螺栓連接板結構的動力學模型可以很好的預示結構的動力學特性，固有頻率的最大相對誤差僅為1.69%，且MAC均大于0.91，滿足工程中對結構動力學建模的需求。

3.1 優化算法對比

本文分別采用遺傳算法[22](Genetic Algorithm,GA)、粒子群算法[23](Particle Swarm Optimization,PSO)和水循環算法(Water Cycle Algorithm, WCA)三種全局智能優化算法對算例中的薄層單元的參數進行修正。為了使修正結果具有可比性，將算法的種群數、粒子群數和雨滴數均設置為30，迭代次數均設置為40次。圖8為分別采用三種優化算法執行模型修正問題的收斂情況對比結果。從圖8中可以看出，相比遺傳算法和粒子群算法，水循環算法的收斂曲線下降最快，目標函數值最小，表明水循環算法在求解結構動力學模型修正問題時具有收斂快、精度高的優點。

圖 8 不同智能優化算法優化結果對比圖Fig.8 Comparison of optimization results of different intelligent optimization algorithms

為了進一步驗證水循環算法在求解最優化問題(式(9))時具有唯一性，本文在保證相同優化參數設定的前提下，重復運行求解過程10次，均得到了相同的最優解和最佳適應度值，即，E1=0.001 MPa；E2=0.2957 MPa；∏(p)=0.006 561 7，說明了采用水循環算法在求解結構動力學優化問題時具有很好的魯棒性。

3.2 抗噪性驗證

本節對所提動力學建模方法的抗噪性進行驗證。通過對實驗模態數據添加人為噪聲來模擬帶有測試噪聲的模態參數識別結果，即：

式中：XE代表固有頻率向量或第i階模態振型向量；η為噪聲水平；δ為均值為0、方差為1，并滿足標準正態分布的隨機向量。

眾所周知，測試噪聲對識別結構的固有頻率和模態振型的影響是不相同的，固有頻率一般可以比較精確地識別，而模態振型的識別精度比固有頻率的精度要差很多。因此，選取了兩種噪聲水平工況來考察本文所提出的螺栓連接板結構動力學建模方法的抗噪性，具體為：工況1(固有頻率( η=1% )，模態振型( η=5%)；工況2(固有頻率( η=2% )，模態振型( η=10%)。表5和表6分別給出了不同噪聲水平下的固有頻率和MAC值對比結果。圖9示出了不同噪聲水平下修正后獲得的MAC值對比情況。

表 5 不同噪聲水平下的固有頻率對比Table 5 Comparison of modal frequencies under different noise levels

表 6 不同噪聲水平下的MAC值對比Table 6 Comparison of MAC under different noise levels

圖 9 不同誤差水平下修正后的MAC值對比Fig.9 Comparison of MAC under different noise levels

綜合上述對比結果可以看出，在第一種噪聲水平下，固有頻率的最大相對誤差僅為2.49%，MAC值均大于0.9，其最小值為0.907。在第二種噪聲水平下，固有頻率的最大相對誤差為4.61%，MAC值均大于0.9，其最小值為0.911。上述結果表明，在兩種噪聲工況下建立的螺栓連接板結構的動力學有限元模型的動力學特性都和真實結構的測試結果吻合得很好，表明本文所建立的動力學建模方法具有良好的抗噪性。

4 結論

本文針對螺栓連接板的動力學建模問題，提出了一種結合薄層單元法和模型修正的建模方法。采用薄層單元法建立了用于模型修正的初始動力學模型，將薄層單元的不確定材料本構參數作為修正變量，引入水循環算法來進行優化求解。從本文的研究中可得出以下結論：

(1)相對于遺傳算法和粒子群算法，采用水循環算法進行動力學模型修正，極大的提高了求解效率。

(2)采用兩端固支五螺栓搭接板結構作為實驗驗證模型，驗證了提出的螺栓連接板結構動力學建模方法是合理可行的，建立的動力學預示模型可以更好地用于進行結構的動力學分析，并且具有良好的抗噪性。