基于非均勻分布復彈簧單元的螺栓連接薄板結構動力學有限元建模

劉曉峰， 孫 偉， 方自文

(1. 東北大學機械 工程與自動化學院， 沈陽 110819；2. 東北大學 航空動力裝備振動及控制教育部重點實驗室， 沈陽 110819)

螺栓連接幾乎應用在所有機械裝備中，其連接界面(結合面)的力學特性對整個機械裝備的靜力學及動力學性能均有重要的影響。因而，在機械裝備設計研發過程中必須重視對螺栓結合部的性能分析與設計。而創建一個簡單、高效并且具有較高精度的結合部模型則是分析和設計螺栓連接結構的重要基礎。

長期以來已有大量學者致力于螺栓結合部的力學建模與分析，真實的螺栓結合部處存在螺栓、螺母等構件，在建模時若將其體現出來，建模過程會較復雜。這類高保真建模方式通常適用于靜力學分析，例如，Al-Nassar等[1]建立了完整的螺栓連接板三維有限元模型并對其在剪切和拉伸載荷作用下的受力性能進行了分析。Gray等[2]提出了一種高效的螺栓連接整體三維模型建立和驗證的方法。為了方便進行螺栓結合部的力學研究，一些建模方法在建模時不體現螺栓、螺母等構件，而是用一些如彈簧、阻尼類的具有集中參數的單元，或是利用一種虛擬介質等來模擬結合面的力學特性。這類簡化模型的建模思路又可分為線性建模和非線性建模。

對于線性建模，通常將模擬結合部力學特性的集中參數或虛擬介質的材料參數等視為常量，例如，Ye等[3]將虛擬材料參數導入有限元分析軟件，建立了包含結合部虛擬材料模擬的螺栓連接結構有限元模型。Zhou等[4]提出了單自由度黏滑摩擦模型，對螺栓連接復合板進行了多尺度建模和有限元分析，并與試驗結果對照，驗證了模型的正確性。McCarthy等[5]用彈簧單元模擬結合部建立了三螺栓連接復合材料梁的分析模型，研究螺栓孔間隙對多螺栓連接結合部載荷分布的影響。

實際螺栓連接結構，由于受激振力變化等因素的作用，其結合部處通常會產生非線性的力學現象，此時線性建模的思想不再適用，需進行非線性建模。如，Zhao等[6]提出了一種用非線性虛擬材料來描述螺栓連接力學特性的方法。Taheri-Behrooz等[7]建立了考慮節點非線性材料特性的單搭接多節點非線性彈簧-質量模型來模擬螺栓結合部。 Kong等[8]提出了一種包含非線性接觸力、阻尼和連接界面幾何特性的螺栓連接結構降階模型，降低了求解非線性動力學方程的計算成本。Armand等[9]提出了一種新的多尺度方法，將螺栓結合部界面粗糙度、接觸壓力和接觸剛度聯系起來，并與多諧波平衡求解器相結合，可以計算不同界面粗糙度下結構的非線性動態響應。

螺栓結合部的非線性建模往往更能再現連接結構真實的振動行為，但是非線性建模的研究通常源于對線性模型的改進。當前，在用彈簧單元模擬螺栓結合部的研究中，大多是用一根彈簧單元或者剛度相等的分布彈簧單元來模擬螺栓及連接區域的剛度。這種建模方法通常無法考慮螺栓影響區剛度非均勻分布對結構的影響。還有其彈簧單元剛度通常是以結構、螺栓材料和螺栓載荷為參數通過經驗計算公式得到，其計算結果的準確性也有待商榷。

螺栓影響區對整個結構系統振動特性的貢獻主要體現為剛度及阻尼。因而，在考慮影響區面積的前提下，本文提出一種用剛度非均勻分布的復彈簧單元來模擬螺栓結合部力學特性的建模方法，并在螺栓連接薄板振動特性分析上進行了實踐。文中詳細描述了上述建模方法的原理，以及利用基于遺傳算法的反推辨識技術確定復彈簧單元參數的流程。最后，以用兩個螺栓連接在一起的薄板結構為例進行了實例研究，通過與試驗的比對證明了上述基于非均勻分布復彈簧單元模擬螺栓結合部力學特性的合理性。

1 基于非均勻分布復彈簧單元的螺栓結合部建模

以下具體描述利用非均勻分布復彈簧單元模擬螺栓結合部力學特性的建模方法，包括建模理念以及復彈簧單元剛度及阻尼參數的確定方法。

1.1 建模理念

螺栓在固定兩個連接結構時會在結合面產生一個以螺栓桿軸心為中心的一個非均勻分布近似圓形的壓力區域，稱為螺栓影響區[10]。在振動分析中，可以說對螺栓結合部的模擬，就是要對這個螺栓影響區的有效模擬，具體涉及確定影響區的面積以及影響區剛度及阻尼特性的分布等。

1.1.1 螺栓影響區的面積

如圖1所示，dw為螺栓承載面直徑，t為連接板的厚度，D為螺栓影響區域直徑。可見，獲得了螺栓影響區直徑就確定了整個影響區的面積。

(a) 螺栓連接板結構

當螺栓型號及連接結構已定時，可以確定dw及t的值，而半頂角α取值范圍一般在25°≤α≤33°[11]，基于上述參數，螺栓影響區直徑的計算式可表達為

D=dw+2t·tanα

(1)

按上述理論計算的螺栓影響區面積為圓形，但在有限元建模時(尤其是在自行研發的有限元程序中)，對圓形區域劃分網格不夠方便，為了簡化建模，這里采用一個與圓形區域面積相接近的正方形區域來描述螺栓影響區的力學特性。因而在本研究中，螺栓影響區的面積僅作為建模時添加復彈簧單元區域的參考，并沒有將實際螺栓影響區面積引入分析模型。

1.1.2 螺栓影響區的剛度及阻尼分布假定

研究表明[12]，在螺栓影響區中靠近螺栓區域的剛度值較大而遠離螺栓區域的剛度值逐漸降低，直至影響區邊界剛度降至0。因而，采用均勻分布的彈簧模擬螺栓影響區的力學特性是不恰當的。本文提出用三種剛度非均勻分布彈簧單元來模擬螺栓影響區的力學特性，包括正弦、拋物線和線性分布，如圖2所示。從圖中可以較為直觀的看出，三種分布均反映了螺栓影響區中心處剛度值最大，中心區向外圍區剛度值逐漸遞減的趨勢，因而可以較為真實地描述結合面的力學行為。

(a) 正弦分布

需要說明的是在本文的建模方法中，彈簧單元僅用于模擬螺栓影響區的面外剛度，且每根彈簧的剛度定義為復剛度k*，該復剛度同時包含了結合面的剛度及阻尼，表達為

k*=k+ic

(2)

這里k表征了螺栓結合面的剛度，而c描述了結合面的阻尼，且這是一種結構阻尼。

首先討論影響區的剛度值分布。基于圖2所描述的螺栓影響區剛度分布假定，在影響區內任一彈簧剛度k可由影響區中心處剛度值kmax以及彈簧單元距離影響區中心距離的2倍(即圖2中的d)來綜合確定。對于正弦、拋物線和線性分布，在影響區內任一彈簧剛度的求解式可分別描述如下

k=kmaxsin[90°(1-d/D)]

(3)

k=kmax[1-(d/D)2]

(4)

k=kmax(1-d/D)

(5)

接著討論螺栓影響區的阻尼值分布。結合面的阻尼分布與剛度分布通常是不一樣的，通常在壓力大的地方剛度大而阻尼小[13]。考慮到阻尼機理的復雜性，為了簡化建模，將每根彈簧的結構阻尼取同一值(即將阻尼在影響區內視為均勻分布)。

綜上，基于非均勻分布復彈簧單元的螺栓結合部建模過程可簡要描述如下：① 根據螺栓的類型確定影響區的面積；② 基于此面積完成螺栓影響區網格劃分，在相應節點上建立同時包含剛度及阻尼的復彈簧單元：③ 指定一種非均勻分布來描述螺栓影響區的彈簧單元剛度，用均勻分布來描述螺栓影響區的阻尼；④ 利用參數辨識技術確定復彈簧單元的剛度及阻尼值。

1.2 復彈簧單元剛度及阻尼參數確定方法

參見式(3)～式(5)，由于假定剛度值滿足一種特定的分布，因而在實際研究中只需要確定中心區最大剛度值kmax，而阻尼值對于每根彈簧單元是一樣，即需要確定參數c的值。

這里采用反推辨識法[14-15]確定上述待定參數，反推法基本思想可概括為：將需辨識的量值定義為設計變量(這里為kmax和c)，通過一個匹配算法不斷修正理論模型中的設計變量值，使計算獲得的結構振動特性參數與試驗測得的相一致，進而反推出待辨識的參數。用于模擬螺栓結合面力學特性的復彈簧單元剛度參數直接影響螺栓連接板結構的固有頻率，而復彈簧單元的阻尼參數直接影響螺栓連接板結構的頻響函數幅值。因而，這里提出通過實測螺栓連接薄板結構的固有頻率和頻響函數幅值來反推復彈簧單元的剛度及阻尼參數。

圖3為利用反推辨識法確定復彈簧單元剛度及阻尼參數的流程，主要包括三項內容，分別是試驗測試、理論計算和匹配計算。對于試驗測試，這里通過錘擊試驗獲得螺栓連接薄板結構的頻響函數，進而由頻響函數峰值處對應的頻率值確定各階固有頻率，這就為反推辨識提供了完整的試驗數據。在理論計算中，這里主要利用自編的有限元程序來完成螺栓連接薄板的有限元建模，詳見第2章。

圖3 復彈簧單元剛度(阻尼)辨識流程

匹配計算的目標是使理論計算結果與實測值盡快匹配，這里采用遺傳算法[16]來實現參數尋優。模型參數的確定是通過將實測的固有頻率(頻響函數值)與相應的有限元模型預測值之間的差異最小化來實現的，故反推辨識彈簧剛度的目標函數可表達為

(6)

(7)

在基于遺傳算法反推復彈簧剛度及阻尼的匹配計算中，首先，需設定好種群數量，即迭代過程中任意一代給出的個體數量(這里個體指復彈簧單元剛度或阻尼參數)、變異概率(種群里個體發生變異行為的概率)、交叉概率(種群里個體發生交叉行為的概率)、迭代次數等參數；接著，實施迭代計算，遺傳算法將通過選擇、交叉和變異等遺傳操作進行概率化的迭代尋優；最后，當達到所設定收斂條件(例如達到最大的迭代次數)后，即可獲得待辨識的復彈簧剛度或阻尼值。

2 螺栓連接薄板結構動力學有限元建模及分析

這里用非均勻分布復彈簧單元模擬螺栓影響區的力學特性，完成螺栓連接薄板結構動力學有限元建模與分析。為了更好地實施所研發的建模理念以及反推辨識確認復彈簧單元的剛度及阻尼參數，這里利用Matlab軟件自行研發有限元程序。圖4為處于懸臂狀態的螺栓連接薄板結構，兩塊薄板由若干個螺栓搭接在一起。以下從板單元的分析、復彈簧單元引入以及螺栓連接板結構振動特性求解描述整個分析的過程。

圖4 螺栓連接薄板結構

2.1 薄板單元分析

選用如圖5所示的4節點薄板單元，每個節點具有3個自由度，即一個位移自由度w(撓度)、兩個轉動自由度θx和θy，單元的長度及寬度分別為a和b。

圖5 4節點薄板單元

單元節點位移矢量為

δe=[w1θx1θy1…w4θx4θy4]T

(8)

每節點的形函數可表示為

(9)

i=1,2,3,4

式中：ξ，η為介于-1和+1間的局部坐標；ξi和ηi為第i個節點的局部坐標。

根據單元剛度矩陣定義可求出單元剛度矩陣，求解式為

Ke=t?BTDBdxdy

(10)

式中：B是由幾何方程推得的單元應變矩陣;D為由材料參數確定的彈性矩陣

相應地，薄板單元質量矩陣求解式可表示為

Me=t?NTρNdxdy

(11)

式中：ρ為板的密度；N為由式(9)描述的節點形函數組成的形函數矩陣。

2.2 復彈簧單元的引入與系統組集

在螺栓影響區的每個節點上都有一個同時包含剛度及阻尼的復彈簧單元，其剛度方向對應于上述板單元的位移自由度為w(撓度)，在有限元分析時，需要將這些單元與板單元進行有效的組集。

復彈簧單元的剛度及阻尼矩陣可分別表達為

(12)

(13)

則總的復彈簧單元復剛度矩陣表達為

(14)

引入邊界條件后，將彈簧的復剛度矩陣與兩塊板的剛度矩陣進行組集，整個組建過程可用式(15)進行描述。式(15)中，左上角框圈起來的部分為第一塊板的剛度組集，m表示第一塊板的節點自由度數，左下角為第二塊板的剛度組集，n表示第二塊板的節點自由度數，中間部分為搭接部分與復彈簧單元復剛度矩陣的組集，l表示復彈簧單元個數。

(15)

通過這種剛度組集方式，得到螺栓連接結構的復剛度矩陣K*，該復剛度矩陣同時包含了結構系統的剛度及螺栓結合面產生的結構阻尼。螺栓連接薄板結構的總質量矩陣M不涉及復彈簧單元，因而總質量矩陣可表達為

(16)

2.3 螺栓連接薄板結構振動特性求解

懸臂狀態下的螺栓連接薄板結構的頻域運動方程為

(K*-ω2M)X=F

(17)

式中：ω為激振頻率;X為響應向量;F為激振力向量。

用于求解螺栓連接薄板結構模態振型的特征方程可表達為

(18)

式中：KR為復剛矩陣K*的實部，ωr和φr分別為第r階固有頻率和實模態振型。

參見式(17)，螺栓連接薄板結構復頻響函數矩陣為

(19)

本文中螺栓連接板結構模型可認為是小阻尼系統，可認為引入的阻尼并不影響系統的特征向量，因而可用實模態理論分析系統振動行為。

利用求得的模態振型矩陣φ，對上述方程做正交變換有

(20)

式中，ηr為由復彈簧單元的阻尼部分產生的第r階模態損耗因子。

在實測時，通常無法也沒有必要獲得整個系統的頻響函數矩陣。只需獲得在i點激勵和j點拾振獲得的頻響函數，這對應頻響函數矩陣的一個元素，可表示為

(21)

實際研究中一般只需要考慮前若干階次即可，即考慮前n0個階次。

3 實例研究

3.1 問題描述

這部分以螺栓連接鋼板為例描述本文提出的用非均勻分布復彈簧單元模擬螺栓結合部力學特性進而完成組合結構動力學建模與分析的方法。圖6(a)為所研究的螺栓連接薄板結構，其由兩塊鋼制帶孔薄板經螺栓連接緊固構成，相關的幾何及材料參數如表1所示。采用錘擊法對該螺栓連接薄板結構進行固有特性測試，對螺栓施加3 N·m的預緊力，測試中通過PCB SN 30272力錘對薄板連接結構施加寬頻激勵，使用Polytec PDV-100激光多普勒測振儀拾振，LMS SCSDAS數據采集分析儀用于獲取激勵及響應信號，在獲取模態振型時，錘擊點為圖6(a)中標注的所有節點，拾振點位置如圖不變，利用獲取的頻響函數可得到各階固有頻率。而文中分析所選取的頻響函數在試驗獲取時的錘擊點位置在圖6(a)中做了突出標注，相關測試結果列在后續的與理論分析相對照的各表格中。

(a) 螺栓連接薄板

表1 螺栓連接薄板結構中板的相關材料及幾何參數

3.2 螺栓連接薄板動力學建模

用于連接兩塊薄板的螺栓為M6外六角頭螺栓，因而參照式(1)可確定每個螺栓影響區的面積為171.37 mm2。所創建模的有限元模型見圖6(b)，每個板共劃分了400個單元，882個節點，兩塊板在螺栓影響區對應節點處用復彈簧單元進行連接。圖6(b)中的節點包含了圖6(a)中的所有節點，以保證分析獲得的頻響函數能與試驗結果比較。考慮到圓形影響區在有限元建模時難以模擬，因而這里用面積相近的矩形區替換實際的圓形影響區來對結合部進行有限元建模。每個影響區共設置了16個板單元，其總的面積與前述圓形影響區面積計算值大體一致。用25個復彈簧單元模擬螺栓結合面的剛度及阻尼特性，且其中復彈簧單元的剛度值共劃分為6個梯度，分別按正弦、拋物線和線性分布賦值。參照式(3)～式(5)，每個復彈簧單元的剛度值都可以通過螺栓影響區中心處最大剛度值kmax計算獲得。

為了明確用哪種復彈簧單元分布形式模擬螺栓影響區力學特性更能提升分析模型的精度，這里進行了對比研究。同時，為了說明這種用分布彈簧單元模擬螺栓影響區建模理念的先進性，還對比了用單根以及均布彈簧單元模擬螺栓結合部力學特性的結果。分別采用反推辨識技術，利用實測的前5階固有頻率，針對不同的螺栓影響區建模方式，確定彈簧單元的剛度值，其中對前5階分配的權重依次為0.1、0.1、0.3、0.3和0.2，遺傳算法中種群數量、變異概率、交叉概率、迭代次數依次設置為50、0.05、0.9和50，辨識得到的螺栓影響區中心彈簧剛度值結果如表2所示。

表2 反推法辨識獲得的模擬螺栓影響區中心處的彈簧單元剛度

需要說明的是表2中的彈簧剛度值是使對應的螺栓影響區模擬方法達到最高模擬精度的剛度值。用此剛度值獲得的螺栓連接板的前5階固有頻率與試驗值的比對分別如表3和表4所示。

表3 單根彈簧和彈簧剛度值均勻分布時模型固有頻率與試驗固有頻率對比

表4 彈簧剛度值非均勻分布時模型固有頻率與試驗固有頻率對比

從表3與表4的對比可非常直觀地看出用單根彈簧或彈簧剛度值均勻分布模擬螺栓影響區時模型獲得的固有頻率與試驗結果的接近程度明顯差于彈簧剛度值按正弦、拋物線和線性分布時的結果。進一步，利用均方根誤差(RMSE)方法對比上述三種非均勻分布模擬時的各階固有頻率與試驗值的偏差，RMSE的求解式如下

(22)

將表4中的數據代入到式(22)，獲得對應三種分布形式的仿真與試驗前5階固有頻率的RMSE，結果如表5所示。

表5 各非均勻分布形式中仿真與試驗前5階固有頻率的RMSE

通過表5中的RMSE對比可以看出，采用正弦分布時，模型前5階固有頻率與試驗值的接近程度更好。故在后續分析時，對于本文的結構，螺栓影響區均采用剛度值按正弦分布的非均勻分布彈簧單元來模擬。

為了進一步描述清楚該建模方法的可用性，以下進一步討論了螺栓影響區面積對分析結果的影響。將螺栓影響區增加至用25個板單元來模擬，相應的復彈簧單元數量增加到36個，并按正弦分布賦彈簧剛度值，求解該搭接板結構的固有頻率。將獲得的結果與上面用16個板單元模擬影響區面積的結果進行比對，如表6所示。從對比結果可以看出，變動影響區面積會使分析結果發生改變，但變動的幅度并不大。可見只要基于一定的理論確定螺栓影響區面積參考值，在實際建模中基于此參考值大致確定一個影響區面積作為添加復彈簧單元的依據則可獲得具有一定精度的分析結果。

表6 影響區面積變動前后模型固有頻率對比

3.3 振動特性求解

前一部分已描述了固有頻率的求解，這里繼續求解螺栓連接板的模態振型。仿真計算獲得的模態振型與實測模態振型的對比如表7所示，可以看出兩者的振型基本一致。

表7 試驗與仿真前5階振型對照

接下來進行頻響函數的計算，計算之前需獲得復彈簧單元的阻尼。前面已述，這里采用阻尼值均勻分布的復彈簧單元模擬螺栓影響區的阻尼特性。參照圖6所描述的錘擊點及拾振點，利用式(21)進行頻響函數計算。首先，采用反推辨識技術，利用實測包含前5階固有頻率的頻響函數曲線，依次分配的權重為0.4、0.1、0.3、0.1和0.1，遺傳算法中種群數量、變異概率、交叉概率、迭代次數依次設置為50、0.05、0.9和50，不斷迭代模型中彈簧阻尼單元的阻尼值，使頻響函數曲線盡可能接近實測曲線，完成迭代后，螺栓影響區任一彈簧阻尼單元的阻尼值為c=3.56×107N/m。接著用此阻尼值獲得最終的頻響函數曲線并與實測值比對，如圖7所示。從圖中也可看出，仿真與實測的頻響函數也有較好的接近。

圖7 實測與仿真頻響函數對比

綜上，從固有頻率、模態振型以及頻響函數的比對中，可以得出結論，本文所提出的用非均勻分布彈簧單元模擬螺栓影響區進而實施有限元建模可實現較高的仿真計算精度。

4 結 論

本文采用非均勻分布復彈簧單元模擬螺栓結合部影響區，以螺栓連接薄板為對象，描述了相關的有限元建模及參數辨識方法，得出如下結論：

(1) 提出用剛度非均勻分布復彈簧單元模擬螺栓影響區壓力非均勻分布的力學特性，并給出了詳細的建模流程及方法。實踐表明用本文假定的正弦、拋物線和線性等非均勻分布復彈簧單元模擬螺栓影響區并建立動力學模型，其建模精度遠高于用單根及均布彈簧單元模擬螺栓影響區的模型。

(2) 提出的復彈簧單元同時包含了剛度及阻尼，這種建模方式既可以較為精確的模擬結合面的力學特性又可簡化建模過程，通過所提出的基于遺傳算法的反推辨識流程可較為精確的確定復彈簧單元的剛度及阻尼參數。

(3) 在有效引入非均勻分布復彈簧單元的基礎上，建立了螺栓連接薄板結構的動力學有限元模型，計算得到的前5階固有頻率與試驗測得的固有頻率偏差均小于4%，求解得到的模態振型也與試驗振型基本一致，仿真和實測的頻響函數也有著很好的接近，從而證明了該建模方法的合理性。需要說明的是，由于模擬螺栓結合部的復彈簧單元的剛度及阻尼是常值，因而本文所創建的動力學有限元模型僅能模擬結構線性振動特性。無法再現由于激振力幅的增加而使螺栓連接結構出現的軟式(或硬式)、諧波失真、跳躍等非線性振動現象。后續，可對復彈簧單元模型加以改進，以適應螺栓連接結構非線性振動研究。