劉昕運(yùn) 吳大林 馬吉?jiǎng)?/p>

(1.陸軍工程大學(xué)火炮工程系， 石家莊 050001； 2.西京大學(xué)機(jī)械工程學(xué)院， 西安 710123)

0 引言

磁流變(MR)阻尼器的內(nèi)部通常填充有特殊的磁流變液體，該液體由非導(dǎo)電液體、微小的鐵磁性粉末和少量添加劑混合而成[1]。這種特殊的流體在磁場(chǎng)作用下表現(xiàn)出粘塑性非牛頓流體的特性[2]。流體的屈服應(yīng)力由磁感應(yīng)強(qiáng)度決定，即為磁流變阻尼器的阻尼力可控的原因[3]。相比于膠質(zhì)阻尼器和其他阻尼器[4-5]，MR阻尼器的滯回環(huán)范圍更大[6]。目前，MR阻尼器發(fā)展迅速，已在車輛工程[7-8]、土木工程[9]、醫(yī)療器械[10]、飛行器[11]和武器工業(yè)[12]等領(lǐng)域被廣泛使用。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，有限元分析(Finite element analysis，F(xiàn)EA)和計(jì)算流體動(dòng)力學(xué)(Computational fluid dynamics，CFD)仿真技術(shù)已成為分析磁流變液體和磁流變阻尼器的主要手段。與傳統(tǒng)的一維解析模型相比，F(xiàn)EA和CFD數(shù)值模擬更加準(zhǔn)確，且能夠獲得更多可視化參數(shù)[13]。文獻(xiàn)[14]在COMSOL中耦合了磁場(chǎng)和流體，并進(jìn)行了實(shí)驗(yàn)和仿真。文獻(xiàn)[15-16]在ANSYS/CFX中用CLL語言耦合了磁流變阻尼器模型。文獻(xiàn)[17]模擬了磁流變阻尼器的磁場(chǎng)分布，并考慮了溫度的影響。以上研究屬于邊界已知的“顯式”主動(dòng)運(yùn)動(dòng)問題，而本文研究的沖擊運(yùn)動(dòng)屬于邊界條件不可知的“隱式”被動(dòng)運(yùn)動(dòng)問題。目前，此類問題在國內(nèi)外尚未見相關(guān)研究和報(bào)道。

本文提出單向耦合數(shù)值模擬方法。首先在COMSOL中對(duì)磁場(chǎng)進(jìn)行有限元分析，然后將磁場(chǎng)結(jié)果導(dǎo)入Fluent中，結(jié)合Six DOF模型和UDF文件對(duì)沖擊過程進(jìn)行分析，通過實(shí)驗(yàn)進(jìn)行驗(yàn)證，以期準(zhǔn)確獲取磁流變阻尼器的力學(xué)特性。

1 數(shù)值模型

1.1 結(jié)構(gòu)原理

圖1 磁流變阻尼器CAD設(shè)計(jì)圖Fig.1 CAD design of MR damper

近年來，磁流變阻尼器出現(xiàn)了多種結(jié)構(gòu)，單出桿式[18]和雙出桿式[19]，有氣體腔[20]和無氣體腔[21]，特殊式[22]等。本文涉及的MR阻尼器用于緩沖艦炮系統(tǒng)中的高速運(yùn)動(dòng)的實(shí)驗(yàn)彈，CAD設(shè)計(jì)如圖1所示。由于外部幾何限制，對(duì)阻尼器結(jié)構(gòu)進(jìn)行了特殊設(shè)計(jì)，如圖2所示，主要由活塞桿、磁流體活塞、空氣彈簧活塞、缸體、限位塊、線圈、壓縮空氣、磁流變液組成?；钊麠U兩端直徑相等，彈性回復(fù)力主要來自空氣壓縮?；钊欣@有兩組線圈，其被外殼密封。導(dǎo)線從活塞桿內(nèi)部通過并連接外部直流電源。當(dāng)線圈通電時(shí)，將產(chǎn)生大致如圖2所示的磁場(chǎng)分布。磁場(chǎng)導(dǎo)致磁流變液中的鐵磁固體顆粒按磁感線形成鏈狀，導(dǎo)致環(huán)形間隙中的部分流體表現(xiàn)出極大的屈服應(yīng)力，流體通過間隙時(shí)會(huì)產(chǎn)生較大的阻尼力。阻尼器主要尺寸如表1所示。

1.2 靜磁場(chǎng)FE模型

有限元模型在多物理場(chǎng)仿真軟件COMSOL/Multiphysics中建立。依據(jù)麥克斯韋方程組[10]計(jì)算電磁場(chǎng)，包括恒定電荷密度的連續(xù)性方程、安培定律和法拉第定律，即

(1)

其中

B=μoμrH

(2)

式中H——磁場(chǎng)強(qiáng)度

J——電流密度

A——磁矢量勢(shì)

B——磁通密度

μr、μo——材料磁導(dǎo)率、真空磁導(dǎo)率

根據(jù)阻尼器的對(duì)稱性，生成二維軸對(duì)稱網(wǎng)格，如圖3所示。計(jì)算域總共有70 375個(gè)四邊形網(wǎng)格單元。為了評(píng)估網(wǎng)格的相關(guān)性，生成另一組更加細(xì)密的網(wǎng)格，共281 500個(gè)四邊形單元。比較兩組網(wǎng)格的磁通強(qiáng)度計(jì)算結(jié)果，輸出曲線幾乎相同。這說明生成的網(wǎng)格已經(jīng)非常精細(xì)，結(jié)果誤差與網(wǎng)格尺寸基本無關(guān)。

圖3 磁流變阻尼器的有限元網(wǎng)格Fig.3 FE mesh of MR damper

采用LORD公司的MRF-140CG型磁流變液，圖4為該磁流變液的特性曲線?；钊^和缸體采用碳素結(jié)構(gòu)鋼1045，活塞桿采用低磁導(dǎo)率金屬。線圈導(dǎo)線采用美國線規(guī)AWG 24標(biāo)準(zhǔn)銅線。部分材料參數(shù)如表2所示。

圖4 MRF-140CG型磁流變液的特性曲線Fig.4 Characteristics of MRF-140CG MR fluid

1.3 磁流變液CFD模型

MR阻尼器的CFD分析需考慮湍流和能量傳遞，其中流體控制方程包括連續(xù)性方程、Navier-Stokes方程和能量方程[23]

表2 材料屬性參數(shù)Tab.2 Material property parameters

(3)

式中ρ——流體密度V——流體速度

fb——體積力p——差壓

μ——?jiǎng)恿φ扯萮——比焓

λ——熱導(dǎo)率Φ——粘性耗散項(xiàng)

Sh——源項(xiàng)T——溫度

1.3.1流體本構(gòu)模型

根據(jù)圖4中磁流變液材料屬性定義磁通密度B和屈服應(yīng)力τy的關(guān)系曲線。使用Hill sigmoidal函數(shù)擬合為

(4)

式中，ε=117 024.8，β=0.967 66，k=1.693 63。

目前，針對(duì)粘塑性流體提出了多種本構(gòu)模型[24-26]。本文采用由雙曲正切函數(shù)定義的流變粘度模型[13-14]，即

(5)

式中α——粘度增長(zhǎng)控制因子，取0.03

ζ——粘度增長(zhǎng)控制因子，取0.1

對(duì)速度沖擊過程進(jìn)行CFD分析時(shí)，必須考慮MR流體的可壓縮性。因?yàn)樽枘崞髟跊_擊瞬間壓力變化非常大，如果不考慮流體的壓縮性將導(dǎo)致壓力計(jì)算出現(xiàn)異常。液體壓力用Tait狀態(tài)方程描述，即

(6)

式中m——密度指數(shù)ρ0——參考密度

E0——參考體積模量

1.3.2空氣彈簧

磁流變阻尼器有空氣彈簧結(jié)構(gòu)，空氣腔存在一定初始?jí)毫Γ渲饕饔檬翘峁┗貜?fù)力，使被壓縮的阻尼器自動(dòng)復(fù)位?？諝鈴椈闪s一般依據(jù)氣體多變過程方程來描述，即

(7)

式中p0——初始?jí)毫?，?.5 MPa

Aair——活塞面積，取1 756.15 mm2

W0——初始體積，取52 684.5 mm3

n——?dú)怏w多變指數(shù)，取1.3

x——活塞位移

1.3.3流體網(wǎng)格定義

流體計(jì)算域仍然采用二維軸對(duì)稱四邊形網(wǎng)格。動(dòng)網(wǎng)格使用Layering層鋪技術(shù)，網(wǎng)格被壓縮至一定程度時(shí)，貼近壁面的網(wǎng)格會(huì)坍縮與合并。而網(wǎng)格被拉伸至一定程度時(shí)，貼近壁面的網(wǎng)格會(huì)分裂出新層。為了順利實(shí)現(xiàn)網(wǎng)格層鋪法，在生成網(wǎng)格時(shí)需考慮Mesh Interface。如圖5所示，活塞兩側(cè)的壁面被定義為剛體移動(dòng)壁面。僅生成流體域網(wǎng)格，間隙處的網(wǎng)格被加密，總共包含了95 403個(gè)四邊形單元。

圖5 流體網(wǎng)格定義示意圖Fig.5 Schematic of fluid mesh definition

1.3.4用戶自定義函數(shù)

磁流變阻尼器的活塞移動(dòng)時(shí)，磁場(chǎng)的影響區(qū)域會(huì)跟著變化。所以，對(duì)于非牛頓區(qū)域的跟蹤和捕捉是必要的。與顯式邊界模型不同，隱式邊界模型的“磁場(chǎng)作用區(qū)域”不能直接定義，其需要依賴移動(dòng)壁面的實(shí)時(shí)坐標(biāo)。因此，基于C語言編寫了用戶自定義函數(shù)(UDF文件)，并在求解器中進(jìn)行編譯。

使用“DEFINE_PROPERTY”宏定義計(jì)算域的動(dòng)力粘度屬性。首先，使用“Get_Domain”宏和“Lookup_Thread”宏獲取移動(dòng)壁面的指針，在此之上用“F_CENTROID”宏訪問移動(dòng)壁面的坐標(biāo)信息并存入數(shù)組；然后寫入式(4)、(5)和磁場(chǎng)FE分析得到的磁通密度。以上函數(shù)中涉及2個(gè)變量，其中坐標(biāo)信息從數(shù)組中提取，剪切速率由“C_STRAIN_ RATE_MAG”宏訪問得到；最后用“C_ CENTROID”宏訪問計(jì)算域的所有網(wǎng)格單元，并把非牛頓粘度函數(shù)賦值于受磁場(chǎng)影響的單元。

活塞只在x方向移動(dòng)，而且存在空氣彈簧力和限位塊，故還需要定義6自由度模塊的UDF文件。根據(jù)頭文件six_dof.h，使用“DEFINE_SDOF_PROPERTIES”宏定義剛體屬性，代碼可以與屬性宏編寫到一個(gè)文件中。使用“SDOF_MASS”宏定義剛體的質(zhì)量，其值取阻尼器負(fù)載的質(zhì)量；使用“SDOFO_1DOF_T_P”宏定義自由度為1個(gè)平移自由度；使用“SDOFO_LOC， SDOFO_MIN，SDOFO_MAX”宏定義剛體的位移限制；最后用“SDOF_LOAD_F_X”宏寫入空氣彈簧力式(7)。

2 數(shù)值結(jié)果

對(duì)于單向耦合模型，式(4)、(5)是耦合磁場(chǎng)FEA和CFD的關(guān)鍵。通過靜磁場(chǎng)分析得到磁通密度B在流體域的分布，然后將其轉(zhuǎn)換為動(dòng)力粘度在流體域的分布，并導(dǎo)入CFD模型。

2.1 靜磁場(chǎng)FE分析

當(dāng)輸入電流I=4.0 A時(shí)，計(jì)算域內(nèi)的磁通密度等值線如圖6所示。大部分磁感線從活塞的兩端穿過流體進(jìn)入阻尼器的缸體，由于流體和活塞桿的磁導(dǎo)率都非常低，磁通量更多地分布在磁性材料中。輸出環(huán)形間隙中間線上的磁通密度分布曲線如圖7所示，磁通密度存在兩個(gè)峰值，分別在環(huán)形間隙的兩端，最大值達(dá)到1.96 T。而在兩個(gè)峰值的中部以及外側(cè)，磁通密度迅速降低至接近0。

圖6 流體域中的磁通密度分布Fig.6 Magnetic flux density distribution in computational domain

圖7 環(huán)形間隙中線上的磁通密度分布Fig.7 Magnetic flux density distribution on middle line of annular gap

2.2 模型驗(yàn)證

模型采用瞬態(tài)壓力基求解器(PBS)。使用Realizablek-ε模型，激活能量方程和粘性熱。設(shè)置Coupled算法并采用二階離散格式，求解控制Courant數(shù)為50，最大求解時(shí)間步長(zhǎng)為1×10-5s。

磁通強(qiáng)度在流體中的主要作用區(qū)域如圖8所示。受磁場(chǎng)影響的主要區(qū)域不全在環(huán)形間隙中，間隙的兩端存在溢出?？赡軙?huì)影響阻尼器的力學(xué)特性計(jì)算結(jié)果。因此，在UDF代碼中編寫環(huán)形間隙兩端的“圓弧”區(qū)域，其結(jié)合圖7的磁通密度分布曲線能夠更準(zhǔn)確地描述非牛頓粘度。

圖8 環(huán)形間隙兩端的圓函數(shù)區(qū)Fig.8 Circular function area at both ends of annular gap

在CSS-55100型萬能試驗(yàn)機(jī)上對(duì)阻尼器進(jìn)行勻速準(zhǔn)靜態(tài)運(yùn)動(dòng)實(shí)驗(yàn)，如圖9所示。固定缸體，對(duì)活塞桿進(jìn)行勻速的加載和卸載。利用壓力傳感器和位移傳感器采集阻力和位移。圖10、11分別為勻速準(zhǔn)靜態(tài)運(yùn)動(dòng)的阻抗力-時(shí)間曲線和滯回曲線。對(duì)比阻尼器的實(shí)驗(yàn)和仿真曲線，驗(yàn)證了模型的準(zhǔn)確性。

圖9 勻速準(zhǔn)靜態(tài)實(shí)驗(yàn)設(shè)備Fig.9 Uniform quasi-static experimental equipment

圖10 勻速準(zhǔn)靜態(tài)運(yùn)動(dòng)的阻抗力-時(shí)間曲線Fig.10 Resistance-time curves of uniform quasi-static motion

圖11 勻速準(zhǔn)靜態(tài)運(yùn)動(dòng)的滯回曲線Fig.11 Hysteresis loop of uniform quasi-static motion

2.3 沖擊CFD分析

當(dāng)剛體初速度為7.4 m/s，壁面邊界的運(yùn)動(dòng)狀態(tài)由流固雙向耦合的6自由度模塊自動(dòng)解算，總仿真時(shí)間為80 ms時(shí)，整個(gè)沖擊過程的等值線和流線圖如圖12所示。沖擊初期，流體在環(huán)形間隙中瞬時(shí)最大速度達(dá)到了75 m/s以上，但隨著動(dòng)能的耗散，流體速度會(huì)很快降低，達(dá)到最大位移僅需4 ms。磁流變阻尼器在運(yùn)動(dòng)過程中，流體產(chǎn)生了許多不停變化的渦，且流線極為不規(guī)則。但較大的渦都出現(xiàn)在活塞運(yùn)動(dòng)方向的后面。

阻尼器吸收的動(dòng)能主要轉(zhuǎn)化為流體熱能。由于液體具有粘性，內(nèi)部摩擦將產(chǎn)生熱量，初始溫度為300 K時(shí)，流體溫度分布如圖13所示。高溫液體主要產(chǎn)生于環(huán)形間隙，因?yàn)榇颂幖羟兴俾屎蛣?dòng)力粘度最大。沖擊運(yùn)動(dòng)結(jié)束時(shí)溫度分布并不均勻，活塞一側(cè)的平均溫度暫時(shí)高于另一側(cè)。

輸出的磁流變阻尼器運(yùn)動(dòng)學(xué)特性曲線如圖14所示。MR阻尼器的壓縮過程持續(xù)4 ms，拉伸過程持續(xù)大約70 ms。圖14a顯示，MR阻尼器最終無法完全復(fù)位，在位移為1.3 mm處幾乎停止運(yùn)動(dòng)，這是由Bingham流體性質(zhì)決定的。阻尼器在復(fù)原運(yùn)動(dòng)時(shí)，流體域的剪切速率逐漸降低。由式(5)可知，在剪切速率很低時(shí)，表觀粘度會(huì)逐漸大幅升高，從而導(dǎo)致流體剪切應(yīng)力變大。而隨著位移行程的減小，空氣彈簧力逐漸降低。當(dāng)空氣彈簧力小于流體阻力時(shí)，活塞開始減速直到停止運(yùn)動(dòng)。阻尼器在復(fù)原運(yùn)動(dòng)階段先加速后減速的特征如圖14b所示。在實(shí)際應(yīng)用中，這種現(xiàn)象并不會(huì)影響磁流體阻尼器的使用。

圖12 磁流變阻尼器速度等值線和流線Fig.12 Velocity contours and streamtraces of MR damper

圖13 磁流變阻尼器溫度等值線Fig.13 Temperature contours of MR damper

圖14 沖擊運(yùn)動(dòng)的運(yùn)動(dòng)學(xué)特性曲線Fig.14 Kinematic characteristic of impact motion

這是由于當(dāng)控制電源斷開后，磁流變液會(huì)恢復(fù)牛頓流體特性，阻尼器將很快完全復(fù)位。圖15為沖擊運(yùn)動(dòng)的滯回曲線，阻尼器的阻抗力曲線變化非常不規(guī)則，最大阻抗力出現(xiàn)在沖擊初期，在行程2 mm處達(dá)到最大阻抗力110 kN，而后迅速降低并逐漸趨于平穩(wěn)。復(fù)原階段中的阻抗力保持在200～400 N極低的水平。該磁流體阻尼器在沖擊過程中幾乎吸收了所有的初始動(dòng)能，其緩沖效能較為優(yōu)異。

圖15 沖擊運(yùn)動(dòng)的滯回曲線Fig.15 Hysteresis loop of impact motion

3 結(jié)論

(1)針對(duì)某磁流變阻尼器提出了結(jié)合磁路有限元分析和計(jì)算流體動(dòng)力學(xué)分析的單向耦合仿真研究方法。在實(shí)驗(yàn)驗(yàn)證勻速準(zhǔn)靜態(tài)運(yùn)動(dòng)模型的基礎(chǔ)上，對(duì)模型進(jìn)行了特定沖擊載荷的流固雙向耦合模擬，獲得了其力學(xué)特性規(guī)律。

(2)在隱式邊界條件的CFD分析中，解決了粘塑性流體本構(gòu)定義、非牛頓粘度區(qū)的跟蹤和捕獲、以及被動(dòng)運(yùn)動(dòng)的動(dòng)網(wǎng)格處理等難點(diǎn)。在磁路FE分析得到雙峰分布的磁通密度的基礎(chǔ)上，用雙曲正切函數(shù)定義了粘塑性磁流變液的非牛頓粘度，并用“圓函數(shù)”準(zhǔn)確捕獲非牛頓粘度區(qū)。模型也考慮了液體的可壓縮性和多變方程描述的空氣彈簧力。

(3)在沖擊運(yùn)動(dòng)過程中，流體域產(chǎn)生大小不一的多個(gè)渦，其中最明顯的渦產(chǎn)生于活塞運(yùn)動(dòng)方向的后側(cè)。磁流變液中的熱量主要產(chǎn)生于環(huán)形間隙，最終在活塞一側(cè)的平均溫度高于另一側(cè)。在持續(xù)供電的情況下，阻尼器將無法完全復(fù)位，但電路控制可以解決這一問題。