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

李 林，曾祥坤

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

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

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

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

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

液 －固耦合（Fluid-Structure Interaction，F(xiàn)SI）問題中，液體的作用力引起結(jié)構(gòu)的變形，同時(shí)結(jié)構(gòu)的大位移又影響流場的形態(tài)，這種類型叫做雙向耦合的液－固耦合分析。雙向耦合的液－固耦合方程是非線性方程，因此需要用迭代的方法得到某一時(shí)刻的解。也就是說，計(jì)算得到的是液－固耦合問題的迭代解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 徑向型液壓襯套動(dòng)特性分析的集總參數(shù)模型

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

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

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

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

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

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

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

橡膠主簧的動(dòng)剛度可以由其靜剛度Ks來表示：

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

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

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

計(jì)算得到的橡膠主簧徑向靜剛度為320 N／mm，橡膠主簧在徑向方向力－位移曲線的計(jì)算結(jié)果與前期研究實(shí)驗(yàn)結(jié)果的對(duì)比見圖3。由圖可見，隨著位移的增加，橡膠主簧靜剛度有所增大；計(jì)算值與實(shí)驗(yàn)值較吻合。取修正系數(shù)為f為1．6，計(jì)算得到橡膠主簧的動(dòng)剛度為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ù)表達(dá)式：

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

3．3 液室的體積剛度

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

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

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

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

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

3．4 橡膠等效活塞面積

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

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

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

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

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

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

將式（5）進(jìn)行拉氏變換，其頻域表達(dá)式為：

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

［1］上官文斌，徐馳．汽車懸架控制臂液壓襯套動(dòng)態(tài)特性實(shí)測與計(jì)算分析［J］．振動(dòng)與沖擊，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］李林．液壓襯套動(dòng)態(tài)特性實(shí)測分析及其集總參數(shù)模型建模研究［J］．振動(dòng)與沖擊，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．