液壓襯套集總參數(shù)模型動態(tài)特性液－固耦合有限元分析

李 林，曾祥坤

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

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

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

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

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

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

應(yīng)力的標(biāo)準(zhǔn)是

位移的標(biāo)準(zhǔn)是

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

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

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

根據(jù)液體的連續(xù)方程和動量方程以及激勵與動反力的關(guān)系，得到液壓襯套的復(fù)剛度為：

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

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

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

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

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

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

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

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

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

計算得到的橡膠主簧徑向靜剛度為320 N／mm，橡膠主簧在徑向方向力－位移曲線的計算結(jié)果與前期研究實驗結(jié)果的對比見圖3。由圖可見，隨著位移的增加，橡膠主簧靜剛度有所增大；計算值與實驗值較吻合。取修正系數(shù)為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 橡膠主簧的徑向阻尼

橡膠主簧阻尼系數(shù)表達式：

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

3．3 液室的體積剛度

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

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

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

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

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

3．4 橡膠等效活塞面積

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

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

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

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

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

液體在慣性通道內(nèi)的動量方程為：

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

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

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

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

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

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

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

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

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

3．6 集總參數(shù)識別結(jié)果分析

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

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

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

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

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

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

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

5 結(jié) 論

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

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

文中提出的研究方法在液壓襯套產(chǎn)品開發(fā)設(shè)計階段具有較高的工程應(yīng)用價值，同樣也適用于其他類型具有液壓阻尼特征元件的設(shè)計優(yōu)化。

［1］上官文斌，徐馳．汽車懸架控制臂液壓襯套動態(tài)特性實測與計算分析［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－h(huán)armonics［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］李林．液壓襯套動態(tài)特性實測分析及其集總參數(shù)模型建模研究［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．