基于精確幾何模型梁單元的螺旋彈簧剛度分析

張 健，齊朝暉，卓英鵬，國樹東

(大連理工大學工業裝備結構分析國家重點實驗室，大連 116024)

螺旋彈簧在機械系統中有著重要影響，尤其在車用發動機配氣機構和懸架系統等復雜機構上起著重要作用[1－2]。螺旋彈簧的剛度特性是設計彈簧的重要評價因素。

在復雜機械系統建模中，將螺旋彈簧視為單自由度振動模型或是采用將彈簧按照圈數分割為集中質量，質量與質量之間通過連接具有當量剛度的無質量彈簧[3]，該方法忽略了彈簧的非線性特性，不能準確地反映彈簧的非線性剛度。采用有限元分析能獲得比較精確的彈簧剛度特性，因而，有限元分析廣泛應用于對彈簧特性的評價。早期建立彈簧有限元模型是利用傳統梁單元建模，將彈簧分割成很多部分，每部分之間通過直梁單元連接，該方法忽略了彈簧本身是空間彎曲結構，故而彈簧模型離散的疏密程度和梁單元在彈簧中的比例都會影響計算的精度和計算速度[4]。文獻[5]應用直梁單元的建模思想建立曲梁單元，證明對于螺旋彈簧結構中使用曲梁單元的精度要明顯高于直梁單元，并且同等精度的情況下使用更少的單元數量。文獻[6]提出了節點自由度為 6的 2節點三維曲梁有限單元模型，在一定程度上是一種高效彈簧建模方法。同時，文獻[7]作為彈簧設計指導手冊，對彈簧設計做出了比較系統的分析，提出了經典理論剛度公式，在其有限元分析案例中采用三維實體單元，這種分析手段方便成熟，適合于對彈簧特性的評價，但并未考慮彈簧彎曲變形對彈簧剛度的影響。文獻[8]中的彈簧模型考慮了彎曲變形，但所建立的曲梁單元只適合于小變形情況。

螺旋彈簧模型應選擇一種幾何非線性分析方法。對于細長結構，近幾年發展起來的幾何非線性分析方法有絕對節點坐標法和“精確幾何模型梁單元”。

絕對節點坐標法是放棄了傳統梁理論中的剛性截面假設，將梁看成一種特殊實體，可精確描述單元剛體位移的變形場，它是由 Shabana[9]首次提出，由于該方法避免了對截面有限轉動插值引起的各種困難，是柔體動力學研究的一個重要進展，受到了國內外學者的廣泛關注[10－11]。

“精確幾何模型梁單元”則是一種摒棄了小位移小轉動假設的幾何非線性分析方法。選取描述單元端部形心位置的位移變量和描述單元端部橫截面有限轉動的角度變量作為獨立的描述變量，在單元域內進行插值離散。而傳統“精確幾何模型梁單元”只有剛性截面假設，分析細長結構的時候會出現“剪切閉鎖”，結合 Bernoulli梁理論和細長梁的幾何關系，可避免對截面轉動矢量進行獨立插值，從而避免傳統方法當中的很多問題[12－15]。

本文基于“精確幾何模型梁單元”提出了一種新的螺旋彈簧剛度分析方法。結合螺旋彈簧的設計構型特點和變形規律，尋找螺旋彈簧的描述參數；基于傳統“精確幾何模型梁單元”幾何非線性分析的方法，同時增加形心線與剛截面垂直假設條件，建立彈簧系統虛功率方程，并通過數值算例進行剛度特性分析，對比與經典理論剛度計算結果、傳統有限元計算結果的差異。

1 螺旋彈簧系統描述

對于任意形狀的圓柱螺旋彈簧，其基本參數包含螺旋線圓柱半徑ρ、高度z和極角θ，以描述彈簧截面形心位置。首先，以右手螺旋彈簧為例，在彈簧底圈中心處，建立總體基為如圖1所示。

圖1 總體基和截面坐標系Fig.1 Global coordinate system and section coordinate system

螺旋彈簧螺旋線的形心矢徑為：

考慮圓柱螺旋彈簧只受到軸向載荷作用時，彈簧產生軸向變形，在變形過程中形心線的位置仍保持螺旋形狀[7]，只是基本參數螺旋線圓柱半徑ρ、高度z和極角θ會發生變化，故高度z、半徑ρ和極角θ是彈簧變形的狀態變量。

對于已有彈簧，往往已知彈簧自由狀態下高度和半徑隨極角變化的構型曲線。只需要通過選擇適當的插值方法將構型曲線離散，便可得到螺旋線上任意一點的形心矢量。

以彈簧的每一圈為一小單元，采用Hermite插值方法，螺旋線圓柱半徑插值函數為：

同時，保證2個樣條單元之間節點處的二階導數連續性。依據這個連續性條件，單元節點處的半徑對初始極角的一階導數表示為：

式中：ρi＝[ρ0ρ1...ρn]；S0和S1均為常數矩陣。

從而可將式(2)改寫為：

式中：

同理，彈簧螺旋線高度為：

彈簧螺旋線極角為：

假設簧絲截面為剛性截面，彈簧在壓縮或者拉伸過程中，不僅彈簧螺旋線的形心矢徑發生改變，而且簧絲截面也會發生相應的扭轉。除形心矢徑的描述變量外，還需增加一個簧絲截面的扭轉角。

插值方法同上，有：

故，螺旋彈簧的完備描述參數為：

2 簧絲截面坐標系的建立

基于剛性截面假設，簧絲截面所做的運動為剛體運動，因而可以用固結于橫截面的截面坐標系描述剛體截面的方位，如圖1所示。截面坐標系的原點在彈簧形心線上，es軸是形心線的切線方向，另外et、eb兩根軸分別為彈簧截面的形心主軸。

螺旋線形心線切矢量為：

卡爾丹描述的截面法向矢量為：

考慮彈簧細長結構特點，采用 Bernoulli梁理論，變形過程中梁截面始終保持與形心線切向垂直，截面法向矢量與螺旋線的形心線切線方向一致。通過式(10)和式(11)，求得卡爾丹角為：

由式(12)和式(13)可知，α、β角完全由彈簧形心線的形狀所決定。因而只需要用第3次定軸轉動的卡爾丹角就可以描述簧絲截面的扭轉。

簧絲截面的其他2個矢量為截面形心主軸，分別為：

3 彈簧單元的慣性力虛功率

剛性簧絲截面的慣性力虛功率包含平動慣性力虛功率和轉動慣性力虛功率兩部分，在彈簧單元域(0，2nπ)內積分可得到彈簧單元慣性力虛功率。

彈簧單元平動慣性力虛功率為：

角速度具有明確的物理意義[16]，有：

式中，截面角速度ω在連體基中的分解為：

考慮到方便對比推導曲率矢量，故將截面角速度在截面坐標系下表達。

采用卡爾丹描述，在總體基下的3個轉軸為：

通過式(11)、式(14)和式(15)將其轉化到截面坐標系下，則3個轉軸為：

根據相繼定軸轉動，彈簧截面的角速度為：

角速度矢量在截面坐標系中的分量為：

彈簧單元轉動慣性力虛功率為：

式中，J為剛截面的慣性矩張量。

4 螺旋彈簧的曲率矢量與軸向應變

螺旋彈簧是空間曲梁結構，可以用曲率矢量表征形心線的彎曲扭轉程度，與角速度相同，曲率也具有明確的物理意義，曲率為：

式中：s0為螺旋彈簧形心線初始弧長坐標；曲率矢量κ在連體基中的分解為：

類比角速度的疊加原理，彈簧單元截面的曲率矢量為：

曲率矢量改寫為：

曲率矢量在截面坐標系中的分量為：

曲率矢量之所以須在截面坐標系下表達，這是由曲梁的本構關系所決定。

剛截面作剛體運動，可將曲梁微元體簡化成一條形心線，梁的應變只是變形前弧長坐標s0和時間的函數。變形后弧長也是這兩個坐標的函數，則曲梁的軸向應變為：

5 彈簧的變形虛功率

在彈簧中取一段長為ds0微元體，該微元體是具有初始曲率的曲梁微元體。

微元體可當成一段剛體，剛體的動力學方程為：

式中：m、f分別為曲梁微元體左端面合力矩和合力；

根據虛功率原理，可得微元體的外力虛功率為：

微元體的慣性力虛功率為：

根據虛功率原理，曲梁微元體的變形虛功率為：

將其沿極角積分可以得到彈簧變形虛功率：

根據截面法向與截面形心線方向相同的假設條件，由式(10)和式(31)可得：

式(37)兩端對時間求一階導數，得：

角速度ω對求一階導數：

將式(40)代入式(39)，有：

曲率矢量κ的一階時間導數為：

將式(42)代入式(41)，則有：

將式(38)和式(43)代入式(36)，可得彈簧的變形虛功率為：

構造曲梁單元的本構關系可以描述為：

整理后，彈簧的變形虛功率為：

6 彈簧剛度的計算

6.1 模型降噪方法

系統動力學方程具有很強的剛性，而高頻彈性振動主要源于應力的快速變化，模型降噪的方法是將應力在一個時間區間（t,t+h）內平均[17]，用平均應力替換應力σ，以此來計算柔體的變形虛功率。

平均應力近似表達為：

則變形虛功率為：

將式(47)代入式(48)，變形虛功率為：

與修改前的變形虛功率相比，增加了慣性項和阻尼項，從而降低了系統的固有頻率，可以獲得更高的計算效率。

6.2 系統虛功率方程

彈簧剛度是通過在彈簧軸向方向上施加測試外載荷獲得，在彈簧系統中外力虛功率包含重力虛功率和測試載荷虛功率。

重力虛功率為：

測試載荷虛功率為：

式中：fn為測試載荷；

系統虛功率方程為：

7 數值算例

算例1.圓柱壓縮螺旋彈簧基本參數：中徑30.4 mm，有效圈數 4，簧絲直徑 4.1 mm，節距為11.94 mm，材料彈性模量2.1×105MPa，泊松比0.3，密度7800 kg/m3。

彈簧模型如圖2所示，彈簧的兩端圓并緊并磨平，支撐圈數為1.5。

圖2 彈簧模型/mmFig.2 The model of spring

彈簧模型的曲率如圖3所示。由圖3可知，彈簧曲率不僅在小單元內部連續，而且在每圈樣條單元的節點拼接處也連續，該種彈簧單元具有很好的光順性。

當彈簧支撐磨平面承受承載力的時候，忽略與支撐圈與活動圈的壓緊，將作用載荷等效到有效圈的端部節點。

圖3 彈簧初始曲率Fig.3 Spring curvature in initial state

彈簧單元底部施加固定約束，頂部平面緩慢施加過程壓力載荷，最大載荷為300 N，載荷方程為：

彈簧各圈高度的變化情況如下圖4所示。

頂部施加壓力載荷 270 N~540 N，在彈簧受壓過程中，彈簧剛度變化情況，如圖5所示。

由剛度變化曲線可知，彈簧在壓縮過程中，彈簧的剛度是漸減的，程序計算剛度值與文獻[7]評價一致。

算例2.圓柱拉伸螺旋彈簧基本參數：中徑18.5 mm，有效圈數10，簧絲直徑2.5 mm，節距為3 mm，材料彈性模量2.1×105MPa；泊松比0.3；密度7800 kg/m3。彈簧模型如圖6所示。

圖4 壓縮過程中各圈高度Fig.4 Each spring height in compression process

圖5 圓柱壓縮螺旋彈簧剛度變化曲線Fig.5 Cylindrical compression coil spring stiffness curve

圖6 彈簧模型/mmFig.6 The model of spring

底部固定約束，頂部施加拉力 100 N~500 N。剛度變化如圖7所示。

由圖7剛度變化曲線可知，拉伸彈簧的剛度是漸增型，程序計算剛度值與文獻[7]評價一致。進一步發現，剛度并不是線性增加，而是非線性的。

算例3.圓柱螺旋彈簧基本參數：中徑157 mm；有效圈數 6；簧絲直徑 13 mm；材料彈性模量2.1×105MPa；泊松比0.3。頂部施加壓力1000 N。將剛度仿真結果與理論值比較。

圖7 圓柱拉伸螺旋彈簧剛度變化曲線Fig.7 Cylindrical tensile coil spring stiffness curve

圓柱螺旋彈簧的經典理論剛度公式為：

式中，G為彈簧材料的切變模量。

計算結果如表1所示。

表1 彈簧剛度的仿真結果與理論值比較Table 1 Comparison of spring stiffness simulation results with theoretical values

由表1對比可知，有限元Solid95計算剛度結果與線性經典理論計算結果接近。對比節距為40 mm的 Solid95(小變形、單元數目為 156456和 7932)與 Beam189(大變形、單元數目為 50)，其計算結果基本一致，可以證明該有限元方法的正確性。同時，本文方法與Solid95計算結果基本一致，可以證明本文方法的正確性。

8 結論

本文將彈簧半徑、高度、極角和截面扭轉角作為螺旋彈簧的描述變量，并根據設計曲線進行插值離散。基于改進后的“精確幾何梁單元”建模方法，通過形心曲線切矢量和簧絲截面扭轉角建立截面坐標系。應用高頻彈性振動降噪方法，得到彈簧系統虛功率方程，最后求解得到彈簧剛度。通過對比傳統有限元模型和經典理論剛度公式的計算結果，驗證了該種螺旋彈簧剛度分析方法的合理性和正確性，并為變剛度彈簧設計做好了理論準備。