基于非線性虛擬材料法的栓接結合部動態特性分析

李 鵬, 楊 成, 劉秀宇, 魏 超, 吳齊運

（北京航天控制儀器研究所, 北京 100039）

0 引言

栓接結合部對航天器的影響表現為使箭體結構及箭載設備的整體剛度降低和阻尼增加, 經栓接形式連接的構件及儀器儀表整體剛度中30% ～50%決定于結合部的剛度特性, 阻尼90% 以上來源于結合部。 栓接結合部的力學行為顯著影響到航天器及武器裝備的結構安全和可靠性, 影響搭載儀器的穩定輸出。 在栓接結合部研究領域, 學者們分別從宏觀和微觀角度對栓接結合部動靜態特性進行了大量研究, 涉及航空航天、汽車工業、數控加工等諸多領域［1-5］。 子結構耦合方法［6-7］是辨識栓接結合部剛度、阻尼的主要方法之一, 然而該方法辨識結果易受測試噪聲影響, 且難以建立栓接結合部剛度、阻尼與接觸面積、預緊力等因素之間的函數關系, 因此辨識結果難以被應用于栓接形式的分析與優化。

由于傳統辨識方法的不足, 國內外學者則從微觀角度去研究栓接結合部的接觸機理。 1966 年,Greenwood 等［5］將微觀表面形貌上微凸體的高度分布看作隨機變量, 提出了著名的基于統計學分析的GW 接觸模型。 GW 接觸模型假設: 微凸體的峰頂是球形, 具有相同的曲率半徑, 微凸體高度服從高斯分布, 表面粗糙度各向同性, 微凸體彈性接觸且變形相互獨立。 1991 年, Majumdar 等［4,8-9］基于機械加工表面具有自仿射分形特征提出了包含分形維數D和分形粗糙度參數G的M-B 模型。M-B 模型指出: 結合面承受的負載越小, 塑性變形的比例越大。 其研究獲得了真實接觸面積的表達式, 并估算了從彈性變形轉變為完全塑性變形的臨界粗糙體尺寸。 研究結果表明: 尺寸較小的微凸體產生塑性變形, 而尺寸較大的微凸體則產生彈性變形。 2010 年, Jiang 等［10］利用分形幾何和赫茲接觸原理對三種不同加工方式（銑削、磨削、刮研）下的粗糙表面形貌、微凸體彈塑性形變以及單個微凸體接觸剛度進行了研究, 利用構造函數方法計算出表面輪廓的相關參數, 理論計算得出的接觸剛度與不同載荷下實驗獲得的剛度是一致的, 證明了其模型的合理性。 溫淑花等［11-12］在對粗糙表面分形理論和球體與平面接觸剛度問題研究的基礎上建立了具有尺寸獨立性的結合部法向和切向接觸剛度分形模型, 通過對該模型的仿真分析得出法向和切向剛度與結合部外載荷的非線性關系。 2012 年, 黃開放等［13］通過計算特定螺栓預緊力下虛擬材料模型各項參數, 驗證了采用虛擬材料方法模擬栓接結合部的可行性。 2013 年, 基于赫茲接觸理論及分形幾何理論, 田紅亮等［14］考慮結合部法向和切向特性之間的相互作用, 推導出虛擬材料的彈性模量、泊松比等材料屬性, 但是沒有考慮虛擬材料的非線性特征。 2017 年,De Filippis 等［15］采用子空間算法推導接觸區非線性回復力, 解釋了飛機機翼與油箱間因螺栓連接導致的模態相互作用現象。

基于分形接觸理論和赫茲基礎理論, 本文推導出了非線性虛擬材料的等效彈性模量、泊松比和密度特征參數與接觸應力的函數關系。 采用有限元方法獲得栓接結合部接觸表面應力分布, 建立栓接結合部非線性虛擬材料模型。 設計L 形試件驗證所提出模型的有效性, 通過對比理論與試驗結果可知, 所提出模型能夠準確地反映栓接結合部的動態特性, 驗證了模型的有效性。

1 栓接結合部非線性虛擬材料模型

1.1 傳統虛擬材料建模方法

目前對結合面剛度與阻尼的研究較多, 但這些剛度與阻尼數據難以直接應用于實際工程中,且與試驗相比精度較低。 剛度與阻尼不是材料的固有屬性, 而彈性模量、剪切模量、泊松比、密度是材料的固有屬性。 因此, 虛擬材料法可以用來研究栓接結合部的接觸問題, 由于栓接結構在結合部處可以看作是由很多微凸體相互接觸而成,該方法是把微凸體接觸的部分看作是一層厚度相同的虛擬材料層, 其模型如圖1 所示。 虛擬材料模型具有建模簡單和較傳統統計模型具有貼近結合部接觸機理、準確性、普遍性等優點, 但機加工表面的一些物理參數譬如粗糙度、波紋度、平面度等與虛擬材料相關, 但是虛擬材料方法忽略了這些因素的影響。

圖1 結合部虛擬材料模型示意圖Fig.1 Schematic diagram of virtual material model at bolted joint

模態分析中, 虛擬材料的固有屬性即彈性模量E、剪切模量Gx、泊松比v和密度ρ。 依據文獻［16］,虛擬材料的特性參數中, 彈性模量和剪切模量只發生在微凸體發生彈性變形階段, 可用如下理論公式計算得出

式（1） ～式（5）中,D為分形維數,G為分形粗糙度參數,ψ為分形維數D決定的域擴展因子,aL為微凸體最大接觸面積,ac為微凸體臨界接觸面積,E'為當量彈性模量,E1、E2為接觸體的彈性模量,v1、v2為接觸體的泊松比,v'為當量泊松比,E*為無量綱的虛擬材料彈性模量,為無量綱的虛擬材料剪切模量,ρ1、ρ2為兩接觸材料的密度,l1、l2為兩接觸表面實體的厚度。

式（1）、式（3）和式（4）中的參數計算公式分別為

式（6） ～ 式（10） 中,Ar為真實接觸面積,為硬度和屈服強度相關系數,為材料特性,H、σ0.2分別為較軟材料硬度和屈服強度,G'為當量剪切模量,G1、G2為兩接觸體材料的剪切模量。

1.2 非線性虛擬材料栓接結合部模型

本文提出的非線性虛擬材料接觸模型是在圖1右側所示的線性虛擬材料模型基礎上將微凸體接觸部分替換為非線性虛擬材料, 通過對非線性虛擬材料相關材料屬性的定義推導及計算, 進而構建成栓接結合部非線性虛擬材料模型。

本文將以最大平截面積a'L為自變量, 通過對非線性虛擬材料彈性模量、泊松比、密度及材料應力的計算, 推導計算出虛擬材料固有屬性。 虛擬材料的彈性模量、泊松比、密度表現出非線性特性, 即隨著應力的改變而改變。 彈性模量、泊松比、密度的非線性體現在兩方面: 隨著結合部整體外載荷改變引起的應力改變, 虛擬材料彈性模量、泊松比和密度隨之改變; 當外載荷一定時,虛擬材料內部因為應力分布不均導致材料各個部分彈性模量、泊松比和密度也各不相同。 本文還推導出了非線性虛擬材料的應力-應變曲線, 進一步證實并說明了虛擬材料的非線性特征, 虛擬材料的本構模型可用于有限元軟件分析, 彌補了商用軟件在接觸力學分析上的不足。

1.3 栓接結合部的彈性模量、泊松比和密度

由于機加工表面具有分形的特性, 因此分形理論已經被國內外許多學者應用到栓接結合的研究中, Wang 等［7］提出了可以用W-M 函數來表征機加工表面特性, 粗糙表面微凸體平截面積的統計學分布函數滿足

式（12）中,a'為單個微凸體的平截面積, 且有a'=2a,a為單個微凸體的實際接觸面積,為最大接觸面的平截面積。 區別彈性變形和塑性變形的臨界接觸面積可寫為

兩個粗糙表面栓接在一起, 實際上可以看成是很多不同大小微凸體的接觸, 如圖2 所示。

圖2 兩個球形微凸體接觸轉化為球體與剛性平面接觸模型Fig.2 Schematic diagram of two spherical asperities contact model transform into spherical-rigid base contact model

圖2 中, 左側有兩個半徑為R1和R2且接觸的彈性球體, 而兩個球體的接觸問題又可以等效轉化為彈性球體和剛性平面的接觸問題。 一個半徑為R的彈性球體, 在外力f作用下與一個剛性平面接觸并產生形變, 彈性球體產生的形變量為d,由赫茲彈性接觸理論可知, 當微凸體只發生彈性變形時, 法向力f和接觸面積a分別為

對于單個微凸體, 可得接觸應變ε的函數

式（15）中,d和R分別為單個微凸體的變形和等效曲率半徑。

單個微凸體的應力函數可表示為

已知單個微凸體的應力和應變函數, 依據彈性模量的定義, 單個微凸體的彈性模量可寫為

由分形理論可知, 單個微凸體的變形量和等效曲率半徑分別為

將式（18）和式（19）代入式（17）中, 單個微凸體的彈性模量可重新表示為

加權平均彈性模量和相應的實際接觸面積可獲得虛擬材料的等效彈性模量, 當微凸體發生塑性變形時, 彈性模量為零; 當微凸體只發生彈性變形, 即平截面積滿足a'c＜a'＜a'L時, 可獲得等效彈性模量

式（21）中,A0為名義接觸面積。 將式（12）代入到式（21）中, 結果為

式（23）中,β為剪切力和預緊力之比,fτ為靜摩擦系數。

依據接觸表面剪切剛度、剪切變形和剪切力之間的關系, 加權平均剪切模量和相應的實際接觸面積可獲得虛擬材料的等效剪切模量, 當平截面積a'滿足時, 對式（23）積分, 已知剪切模量和相應的實際接觸面積可獲得虛擬材料的等效剪切模量

栓接結合部由表面微凸體組成, 結合部的質量只是很小的一部分。 因此, 虛擬材料的密度對栓接結合部的動態特性影響較小。 依據真實接觸面積與名義接觸面積之比, 提出一個簡化等效密度定義

1.4 栓接結合部的接觸應力

假設: 在整個名義接觸面積A0上, 接觸應力均勻分布, 則平均接觸應力σV為平均彈性微凸體接觸應力σt與平均塑性微凸體接觸應力σs之和。

當平截面積滿足a'＞a'c時, 即微凸體發生彈性變形時, 此時單個微凸體的法向力可表示為

將式（12）、式（27）代入式（28）中, 可得

將式（12）、式（30）代入式（31）中, 可得接觸表面的塑性應力

最終, 可得非線性虛擬材料的平均接觸應力

2 非線性虛擬材料的驗證

2.1 非線性虛擬材料模型試驗驗證

為了驗證所提出模型的準確性, 考慮航天器常用的栓接形式, 本文設計了L 形組合件, 尺寸參數如圖3 所示, 所用L 形試件厚度為20mm。 設計的L 形組合件具有不對稱的特性, 在模態試驗中通過求解振型及各階頻率, 可以很好地體現結合部特性。 L 形組合件為兩個相同尺寸、相同加工方式的L 形件通過螺栓裝配組成, 如圖4 所示, 圖中藍色部分表示的是根據虛擬材料模型理論添加的可在有限元模型中表示結合部特性的結構層。在結合部垂直方向均勻布置2 個φ16 的孔, 其編號為1 和2, 該L 形件分形維數D=1.28, 分形粗糙度參數G=3.5 ×1012, 域擴展系數ψ=2.04, 復合彈性模量E*=92GPa。

圖3 L 形試件尺寸Fig.3 Dimensional diagram of L-shaped specimen

圖4 L 形試件組合體Fig.4 Schematic diagram of L-shaped specimen assembly

L 形件的材料為球墨鑄鐵, 型號為QT600-3,通過磨削加工獲得粗糙表面, 其粗糙度為Ra=1.6μm, L 形件的材料屬性如表1 所示。 為了精確獲得栓接結合部的動態特性, 試驗過程中L 形組合件被繩索吊起以模擬自由模態, 利用8 個加速度傳感器（型號為PCB Model330B30）分別測試試件三個方向的振動數據, 圖4 展示了單個試件上傳感器的布置情況, 另一個試件的傳感器布置與其對稱。采用力錘提供激振力, 連續敲擊8 次, 將振動數據和力錘數據采集并利用LMS Test. Lab 振動噪聲測試系統進行數據分析, 得到L 形組合件各階模參數和振型。 試驗過程中為最大限度消除實際噪聲等各種外在因素的影響, 提高信噪比, 在一定的預緊力情況下, 試驗反復重復測試6 組, 取加權平均值作為最終的實際數值。

表1 L 形件的材料特性參數Table 1 Material properties of L-shaped specimen

（1）傳統虛擬材料模型試驗驗證

在傳統線性虛擬材料方法中, 虛擬材料層的厚度應保持不變, 本文選擇厚度為0.1mm 的虛擬材料層作為研究對象。 在試驗過程中, 兩個螺栓提供相同的預緊力F=10kN, 采用傳統虛擬材料模型可以計算出虛擬材料的平均應力、等效彈性模量、泊松比和密度, 如表2 所示。

表2 傳統虛擬材料法的特征參數Table 2 Characteristic parameters of traditional virtual material method

利用ANSYS 有限元軟件對L 形組合件進行建模和靜力學分析, 獲得栓接結合部各節點壓強值。通過Matlab 軟件, 根據各節點的壓強值計算出各節點對應接觸剛度和接觸阻尼值, 通過MATRIX27矩陣單元建立栓接結合面節點-節點的接觸剛度和接觸阻尼單元。 結合面之間采用目標面單元為TARGE170 單元、接觸面單元為CONTA174 單元進行柔-柔面接觸, 選擇初始綁定, 其他條件默認。

將有限元仿真分析得到的固有頻率與力錘敲擊試驗所得到的固有頻率進行對比, 如表3 所示。結果表明: 傳統虛擬材料方法與試驗結果的誤差均在10%以上, 誤差較大, 傳統線性虛擬材料模型不能夠準確地表征栓接結合部的動力學特性。該方法的不足之處在于沒有考慮螺栓分布對結合部動態特性的影響, 忽略了栓接結合部表面非線性的特點。

表3 L 形組合件均布預緊力下試驗模態與普通虛擬材料模態Table 3 Test modal and common virtual material modal for L-shaped specimen assembly under uniformly distributed preload

（2）非線性虛擬材料模型試驗驗證

栓接結合部的動態特性影響因素不僅僅包括虛擬材料的特征參數, 還包括虛擬材料的厚度。本文引用等效剛度方法［17］來獲得虛擬材料的厚度,假設在整個名義接觸面積A0上的平均壓強分布為PE, 力F和虛擬材料變形δ及栓接結合部剛度的函數關系可表示為

式（34）中,E為虛擬材料的等效彈性模量,h為虛擬材料的厚度,KN為栓接結合部的法向接觸剛度。 根據文獻［18］, 其表達式為

依據式（34）, 虛擬材料厚度為

虛擬材料的厚度和彈性模量成正比關系, 和接觸剛度成反比關系。 根據等效剛度方法, 可以推出虛擬材料厚度和接觸壓強的函數關系, 如圖5所示。

圖5 虛擬材料厚度和壓強的關系Fig.5 Relationship between virtual material thickness and contact pressure

由圖5 可知, 虛擬材料的厚度隨著接觸壓強的增大而增加。 為了建立栓接結合部非均勻壓強分布的模型, 整個接觸表面的虛擬材料層厚度應該保持一致, 故選擇厚度為0.1mm 的虛擬材料層作為研究對象。

2.2 非線性虛擬材料模型有限元分析

對L 形組合件進行靜力學分析, 獲得的栓接結合部的應力分布情況如圖6 所示。

圖6 結合面的應力云圖和AB 線上的應力分布Fig.6 Stress nephogram and stress distribution on AB line of L-shaped specimen assembly bolted joint

從栓接結合部中選擇一條直線AB, 當每個螺栓預緊力為10kN 時, 栓接結合部的應力分布呈現不均勻的特性, 靠近螺栓孔的位置應力最大為5.1MPa, 而在遠離螺栓孔的A 點接觸應力值最小為3.2MPa。

在理想情況下, 結合部虛擬材料由于具有非線性, 因此材料處處彈性模量均不相等。 如果能利用所得到的連續且各點不相等的彈性模量進行模態分析, 則可完全達到分析非線性虛擬材料模態情形的效果。 但由于受有限元分析軟件的限制,本文設計了折中方案, 即根據集合面靜力分析得出的應力分布情況將中間虛擬層劃分為數個區域。在試驗過程中, 由于每一個螺栓的預緊力是相同的, 因此壓力分布是對稱的。 為了研究栓接結合部虛擬材料的非線性特性, 把栓接結合面的一半分成7 個區域, 如圖7 所示, 以此來近似逼近材料非線性效果。

圖7 結合面的虛擬材料區域劃分Fig.7 Virtual material areas division of L-shaped specimen assembly bolted joint

在有限元模型中, 根據栓接結合部非線性虛擬材料不同區域的壓強值可以獲得對應的非線性虛擬材料參數, 如表4 所示。 最后進行栓接結合部結構動力學仿真分析, 獲得栓接結合部的動態特性, 并且通過與傳統虛擬材料法分析結果進行對比, 驗證模型的準確性。

表4 L 形試件結合部虛擬材料的特征參數Table 4 Characteristic parameters of different L-shaped specimen bolted joint virtual material areas

將有限元仿真軟件分析得到的1 ～3 階振型圖與力錘敲擊試驗所得到的振型圖進行對比, 結果如表5 所示。 可以發現, 仿真結果和試驗結果各階振型均高度吻合。

表5 L 形組合件有限元及試驗振型Table 5 Finite element analysis and test vibration mode of L-shaped specimen assembly

2.3 結果對比

傳統線性虛擬材料法和非線性虛擬材料法的試驗結果和理論結果如表6 所示。 非線性虛擬材料法三階模態誤差分別為3.10%、8.80%、2.20%,其誤差小于傳統線性虛擬材料法12.98%、12.60%、-10.50%的結果, 表明非線性虛擬材料法可以更好地表現出結合部在受外載荷情況下的模態特性。

表6 L 形組合件均布預緊力下試驗模態及普通虛擬材料、非線性虛擬材料模態Table 6 Test modal, common virtual material modal and nonlinear virtual material modal for L-shaped specimen assembly under uniformly distributed preload

3 結論

考慮到栓接結合部的動態特性對航天器結構可靠性以及箭載儀器儀表的輸出精度都有重要的影響, 故本文基于分形理論提出了一種基于非線性虛擬材料法的栓接結合部模型。 通過對模型研究, 可得出以下結論:

1）本文將傳統中對結合部的剛度、阻尼的研究轉化為對結合部的非線性虛擬材料固有屬性的研究, 該模型將結合部看作由結合部上部分、結合部下部分及中間非線性虛擬材料層組成, 基于分形理論和赫茲接觸理論, 計算出非線性虛擬材料的等效彈性模量、泊松比和密度特征參數, 并研究接觸表面應力與虛擬材料特征參數的函數關系。

2）本文利用ANSYS 有限元軟件, 通過栓接結構靜力學分析獲得栓接結合部各節點壓強值。 通過Matlab 軟件, 根據各節點的壓強值計算出各節點對應的接觸剛度和接觸阻尼值, 通過MATRIX27矩陣單元建立栓接結合面節點-節點的接觸剛度和接觸阻尼單元。 在有限元模型中, 根據栓接結合部非線性虛擬材料不同區域的壓強值可以獲得對應的非線性虛擬材料參數, 進行栓接結合部結構動力學仿真分析, 獲得栓接結合部動態特性。

3）設計L 形組合件進行試驗, 驗證非線性虛擬材料的可行性。 與傳統線性虛擬材料法進行對比可知: 在相同螺栓預緊力的條件下, 非線性虛擬材料固有頻率的最大誤差為8.80%, 最小誤差為2.20%, 結果遠遠優于傳統線性虛擬材料法10%以上的誤差, 表明非線性虛擬材料模型可以更準確地表征栓接結合部的動態特性。