液壓襯套集總參數模型動態特性液－固耦合有限元分析

李 林，曾祥坤

（1．廣東輕工職業技術學院 汽車系，廣州 510300；2．華南理工大學 機械與汽車工程學院，廣州 510641；3．廣東技術師范學院 機電學院，廣州 510665）

徑向型液壓襯套是汽車底盤隔振系統中廣泛應用的新型隔振器，極大程度地改善了車輛操縱穩定性和平順性［1－2］。液壓襯套非線性動態性能的研究以集總參數模型分析為主，但動特性的計算結果是否與實驗結果吻合，很大程度上取決于模型中物理參數的準確性［3］。目前，這些物理參數的獲得大多依賴于實驗測試或近似解析公式，這在識別精度及預測成本等方面有較大的局限性［4－5］。因此，提出一種不依賴實驗樣件制作的液壓襯套集總參數及動特性的高效辨識預測方法尤為重要。

隨著計算機及有限元技術的發展，液－固耦合有限元非線性分析方法采用 ALE（Arbitrary Lagrangian Eulerian）描述來處理液體－固體耦合邊界的大位移運動，使直接預測計算液壓襯套非線性動力學特性成為可能［6－11］。近年來，國內外一些學者采用非線性有限元技術主要對動力總成液阻懸置的液壓阻尼特征進行了相關研究。Foumai等［6］采用有限元分析方法得到慣性通道式液阻懸置橡膠主簧的剛度、體積柔度和等效活塞面積，但文中僅給出了橡膠主簧軸對稱結構的有限元模型、橡膠材料的超彈性本構模型和計算得到的橡膠主簧的力－變形曲線。Muller等［7］利用ANSYS和FLUENT有限元商業軟件配合，計算一液壓懸置軸對稱橡膠主簧的體積柔度的方法，但該方法只能進行單向液－固耦合有限元分析，且文中未提及其他參數的識別方法。上官文斌等［12－14］在實驗測試的基礎上，研究了慣性通道－活動解耦盤式液阻懸置及其橡膠主簧靜、動態特性的有限元分析方法，對液－固耦合有限元分析中有關的一系列計算和建模方法進行討論，探討了利用有限元仿真技術辨識集總參數的可行性。目前發表的文獻，還未應用非線性有限元分析技術對液壓襯套集總參數識別及動力學特性預測展開全面研究。

文中以徑向型液壓襯套為研究對象，建立其動力學特性分析的集總參數模型，基于ADINA有限元軟件平臺，采用雙向耦合的液－固耦合有限元分析方法辨識模型中的主要物理參數，并對其動態特性進行計算分析。探討了低成本前提下精確預測液壓襯套集總參數模型的參數及動態特性的可行性。

1 雙向液－固耦合的有限元分析理論

液 －固耦合（Fluid-Structure Interaction，FSI）問題中，液體的作用力引起結構的變形，同時結構的大位移又影響流場的形態，這種類型叫做雙向耦合的液－固耦合分析。雙向耦合的液－固耦合方程是非線性方程，因此需要用迭代的方法得到某一時刻的解。也就是說，計算得到的是液－固耦合問題的迭代解X1、X2、…，根據應力、位移或者兩者的結合來檢查迭代的收斂性［15］。

應力的標準是

位移的標準是

式中和分別表示液體和結構的位移和分別為液體和結構的應力，下劃線表示這些值只定義在液－固耦合界面上，n為耦合邊界上的法線方向；ετ、εd分別是應力和位移收斂的容許誤差，ε0是事先給定的常數（≡10－8）。

2 徑向型液壓襯套動特性分析的集總參數模型

徑向型液壓襯套液室及慣性通道的布置情況見圖1（a）。根據徑向型液壓襯套的工作原理，建立其集總參數模型見圖1（b）。模型中，Kr與 Br分別為橡膠主簧的剛度和阻尼系數；Ap為活塞的等效面積；Ii、Ri分別為慣性通道中液體的慣性系數和慣性通道對其中液體流動的流量阻尼系數。KV＝Kv1＋Kv2，BV＝Bv1＋Bv2，其中Kv1與Kv2（單位為N／m5）表示兩液室的體積剛度，Bv1與Bv2（單位為N·s／m5）為兩液室的體積膨脹阻尼。KVV和BVV分別為 KV和 BV的等效線剛度和等效黏性阻尼［1］。

根據液體的連續方程和動量方程以及激勵與動反力的關系，得到液壓襯套的復剛度為：

液壓襯套動剛度，滯后角 φ＝arctg（K1／K2），阻尼系數 C＝Kdsin（h）／ω。

圖1 徑向型液壓襯套Fig．1 The radial hydro-bushing

3 液壓襯套集總參數模型有限元分析方法的參數識別

集總參數模型中的主要物理參數 Kr、Br、Kv1、Ap、Ii、Ri一般由實驗方法或者由近似解析公式計算得到，但是通過這些方法獲得的參數精度，與公式變化范圍、實驗條件和實驗人員等因素有關，而且實驗花費成本較高、周期較長，并且液壓襯套內部流體的運動規律、壓力分布等情況實驗或解析公式是無法獲取的。為了提高參數識別的精度，降低成本，文中采用非線性液－固耦合有限元分析方法辨識模型中的集總參數。

3．1 橡膠主簧的徑向剛度Kr

橡膠主簧的動剛度可以由其靜剛度Ks來表示：

式中：f為修正系數，取值范圍為1．2～1．6之間。

圖1（a）典型徑向型液壓襯套的結構圖中，橡膠主簧的內外表面分別與金屬內管、外管硫化在一起，起支撐和傳力的作用。與橡膠的大變形相比，金屬的變形可以忽略不計，不能影響主要流場的特性。按照有限元分析建模的簡化原則，有限元模型中只考慮橡膠件對變形影響比較顯著的特征，對變形影響較小的特征進行簡化。金屬外管固定安裝在車架上，因此，在進行邊界條件定義時，令該面上所有節點的位移為零。

圖2為橡膠主簧的簡化模型。橡膠材料選用超彈性本構關系中的Mooney－Rivlin模型，橡膠主簧內表面施加位移載荷。網格剖分采用Delaunay法，選用四面體的4／1單元，單元的最大尺寸為5 mm。橡膠主簧的有限元模型共有18 917個單元、4261個節點。

計算得到的橡膠主簧徑向靜剛度為320 N／mm，橡膠主簧在徑向方向力－位移曲線的計算結果與前期研究實驗結果的對比見圖3。由圖可見，隨著位移的增加，橡膠主簧靜剛度有所增大；計算值與實驗值較吻合。取修正系數為f為1．6，計算得到橡膠主簧的動剛度為512 N／mm。

圖2 橡膠主簧的簡化模型Fig．2 The simplified model of rubber spring

圖3 橡膠主簧徑向力－位移特性曲線Fig．3 Radial force-displacement curve of rubber spring

3．2 橡膠主簧的徑向阻尼

橡膠主簧阻尼系數表達式：

其中：M為橡膠主簧的質量，ξr為橡膠材料的阻尼比，實驗得到硬度為邵氏50度的天然橡膠的阻尼比為0．069，計算得到橡膠主簧的阻尼系數為138 N·s／m。

3．3 液室的體積剛度

在外部激勵作用下，液壓襯套金屬內管相對金屬外管發生位移，引起橡膠主簧膨脹變形時，液室內部壓力變化與體積變化之比稱為液室體積剛度。為計算體積剛度構建的固體模型如圖4（a）所示。固體有限元模型中凹陷內腔表面為與液體相接觸的面，定義為液－固耦合面，外管表面所有節點位移為零。

圖4 計算液室體積剛度的有限元模型Fig．4 Calculated finite element model of chamber volume stiffness

圖4 （b）為計算液室體積剛度所需的液體模型，液體的密度和黏度視為常數，定義為不可壓縮流體。液體與橡膠襯套接觸的內弧面定義為液－固耦合面；液體外表面與液壓襯套外管接觸，形成液室，定義為剛性不可滑移的壁面，同時在此面上加載速度載荷。

液體模型外表面施加均勻的速度載荷后，由于液體是不可壓縮的，因此液室產生壓力，使得橡膠襯套產生膨脹變形。通過液－固耦合有限元分析可以求出液室的均布壓力；由加載速度、液體外表面的面積和加載時間可以求出流進液室的液體體積。對于不同加載速度，可以得到液室的壓力與體積的關系，進而求出液室的體積剛度。計算得到液室的壓力變化和體積變化之間的關系曲線見圖5，液室的體積剛度計算值為57×108N／m5。

圖5 液室的體積與壓力的關系Fig．5 Volume-pressure of the chamber

3．4 橡膠等效活塞面積

外部激勵作用下，橡膠主簧發生徑向振動位移x時，液室內部液體體積的變化量等于等效活塞面積Ap所排出的體積。即：

計算橡膠等效活塞面積的幾何模型與計算液室體積剛度的幾何模型一致，有限元模型的定義也基本相同，只是載荷條件有所區別。固體有限元模型內管施加徑向位移載荷x后，可以求出橡膠內腔的徑向位移xp，由xp和液體外表面的面積即可求出液室體積的變化，即為所排出液體的體積。

記錄不同徑向位移載荷作用下xp的變化，計算得到的橡膠等效活塞面積和橡膠徑向位移之間的關系見圖6。由圖可見，當橡膠的徑向位移大于2 mm后，等效活塞面積基本保持不變，其大小與激勵的振幅和頻率無關。計算得到的等效活塞面積為4．13×10－2m2。

圖6 橡膠等效活塞面積Fig．6 The equivalent rubber piston area

3．5 慣性通道中液體的慣性系數和流量阻尼系數

液體在慣性通道內的動量方程為：

將式（5）進行拉氏變換，其頻域表達式為：

因而只要得到慣性通道兩端壓力差ΔP與流經慣性通道的液體流量之間的幅頻特性，就可以識別出流體流經慣性通道時的慣性系數和流量阻尼系數。

計算所需的固體模型與圖2所示的模型相同，液室內壁定義液－固耦合面。液體模型由兩個液室和慣性通道內的液體組成，在一側液室的外弧面上施加壓力載荷；液室及慣性通道的外表面與外管接觸，定義為剛性不可滑移的壁面；其他的面與橡膠接觸，定義為液－固耦合面，如圖7所示。

對圖7所示的有限元模型進行液－固耦合分析時，液體模型所施加的壓力載荷的大小與時間成正弦曲線變化關系，記錄不同頻率下慣性通道另一端流量的變化，根據式（6）計算可得液體流經慣性通道時的質量慣性系數Ii和流量阻尼系數Ri。

由圖8可見，當慣性通道兩端壓力差的變化頻率較低時，液體流經慣性通道時的慣性系數和流量阻尼系數基本變化不大，識別得到的參數 Ii＝1．56×104kg／m4，Ri＝2．5×106N·s／m5。

圖9為激振頻率40 Hz時慣性通道內液體流速的分布情況。隨著頻率的升高，慣性系數和流量阻尼系數都有增大的趨勢，液體在慣性通道內的流動變得困難。這就導致液壓襯套壓力增加一側的液室內部液體無法自由流向另一側壓力較低的液室，液壓襯套傳遞的力相應增加，動剛度增大。

圖7 液體有限元分析模型Fig．7 The finite element analysis model of the fluid

圖8 不同頻率下慣性通道慣性系數和阻尼系數Fig．8 Inertia and damping coefficients of inertia channel

圖9 慣性通道中液體流速Fig．9 Fluid velocity in the inertia channel

3．6 集總參數識別結果分析

基于最小二乘參數估計的方法［16］和基于液－固耦合有限元集總參數識別的方法得到的結果以及兩者之間的相對誤差分析見表1。從表中可以看出，兩種識別方法得到的大部分計算結果比較吻合，誤差較小；只有橡膠主簧阻尼系數的識別誤差相對較大，這主要是由于有限元計算過程中模型簡化造成的。

表1 兩種參數識別方法結果及誤差分析Tab．1 The results and error analysis of two parameter identification methods

4 徑向型液壓襯套動態特性的有限元計算分析

液壓襯套動態特性的振幅相關性和頻率相關性具有較強的耦合特征。采用圖2所建立的固體黏彈性模型和圖7液體模型對其動特性進行有限元分析，固體模型的分析類型選擇動態顯式算法，液體模型設置為瞬態分析模式。在進行動態特性計算時，不考慮溫度及其對液體黏度、密度的影響，在橡膠襯套內管表面上施加徑向諧波位移激勵，其表達式為：

圖10 徑向型液壓襯套動特性計算結果（振幅0．1 mm）Fig．10 Calculated dynamic characteristics of the radial hydro-bushing

式中：A、f分別為動態位移激勵的振幅和頻率。

采用液－固耦合有限元方法計算得到的徑向型單慣性通道式液壓襯套的動剛度和滯后角見圖10。通過與實驗曲線的對比分析，其計算誤差較小，而且較準確的預測了滯后角峰值出現的頻率。與最小二乘參數估計得到的液壓襯套的動態特性曲線相比較［16］，有限元法的計算結果與實驗結果更加接近，因而可以認為利用有限元方法識別得到的液壓襯套動態特性集總參數具有更高的精度。

5 結 論

（1）介紹了徑向型液壓襯套動特性分析的集總參數模型，建立了相應集總參數辨識的有限元模型，探討了運用雙向耦合的液－固耦合有限元分析方法識別集總參數模型中主要物理參數的計算方法和過程，對比分析了計算值與部分實驗結果及最小二乘參數估計的結果，計算結果誤差較小；

（2） 聯合集總參數模型及液－固耦合有限元分析技術，對液壓襯套動態特性進行仿真計算，通過與實驗值及最小二乘擬合結果的對比，結果表明有限元法能準確的預測液壓襯套的動力學性能，且其計算結果較最小二乘擬合結果具有更高的精度。

文中提出的研究方法在液壓襯套產品開發設計階段具有較高的工程應用價值，同樣也適用于其他類型具有液壓阻尼特征元件的設計優化。

［1］上官文斌，徐馳．汽車懸架控制臂液壓襯套動態特性實測與計算分析［J］．振動與沖擊，2007，26（9）：7－10．SHANGGUAN Wen-bin， XU Chi． Measurement and calculation analysis of dynamic characteristics for hydraulic bushing in controlling arm of vehicle suspension［J］．Journal of Vibration and Shock，2007，26（9）：7－10．

［2］Arzanpour S，Golnaraghi M F．Development of an active compliance chamber to enhance the performance of hydraulic bushings［J］．Journal of Vibration and Acoustics，2010，132：1－7．

［3］Gil-Negrete N．Predicting the dynamic behaviour of hydro-bushings［J］．Shock and Vibration，2002，12：91－107．

［4］Sauer W，Guy Y．Hydro bushings innovative NVH solutions in chassis technology［R］．SAE 2003－01－1475，2003．

［5］Arzanpour S．Active and semi-active bushing design for variable displacement engine［D］．Canada：University of Waterloo，2006．

［6］Foumani M S，Khajepour A，Durali M．Application of shape memory alloys to a new adaptive hydraulic mount［R］．SAE 2002－01－2163，2002．

［7］Muller M，Eckel H G．Reduction of noise and vibration in vehicle by an appropriate engine mount system and active absorbers［R］．SAE 960185．

［8］Christopherson J，Jazar G N．Dynamic behavior comparison of passive hydraulic engine mounts．Part 2：Finite element analysis［J］．Journal of Sound and Vibration，2006，290：1071－1090．

［9］Foumani M S，Khajepour A，Durali M．Application of shape memory alloys to a new adaptive hydraulic mount［R］．SAE 2002－01－2163，2002．

［10］Yoon J Y，Singh R．Indirect measurement of dynamic force transmitted by a nonlinear hydraulic mount under sinusoidal excitation with focus on super－harmonics［J］．Journal of Sound and Vibration，2010，329：5249－5272．

［11］ Garcia M J，Kari T L，Vinolas J，et al．Frequency and amplitude dependence of the axial and radial stiffness of carbon－black filled rubber bushings［J］．Polymer Testing，2007，26：629－638．

［12］Wang L R，Lu Z H，Hagiwara I．Integration of experiment and hydrostatic fluid－structure finite element ananlysis into dynamic characteristic prediction of a hydraulically damped rubber mount［J］．International Journal of Automotive Technology，2010，11（2）：245－255．

［13］Shangguan Wen-bin，Lu Zhen-hua．Experimental study and simulation of a hydraulic engine mount with fully coupled fluid-structure interaction finite element analysis model［J］．Computers and Structures，2004，82：1751－1771．

［14］Zhang Yun-qing，Shangguan Wen-bin．A novel approach for lower frequency performance design of hydraulic engine mounts［J］．Computers and Structures，2006，84：572－584．

［15］ Adina R，Adina D．Theory and Modeling Guide-ADINA-A and-F［M］．2010．

［16］李林．液壓襯套動態特性實測分析及其集總參數模型建模研究［J］．振動與沖擊，2013，32（22）：183－188．LI Lin．Experimental study and lumped parameter modeling analysis of dynamic characteristics for hydraulic bushing［J］．Journal of Vibration and Shock，2013，32（22）：183－188．