祝世興，楊麗昆，魏 戩，祝恒佳

（中國民航大學航空工程學院，天津 300300）

磁流變阻尼器（magnetorheological damper，MRD）是一種阻尼可調的半主動智能控制器件，其傳動介質磁流變液（magnetorheological fluid，MRF）是一種受磁場影響的特殊材料［1］。MRF的黏度可隨外加磁場強度不同而迅速發(fā)生變化，進而改變輸出阻尼力的大小。因此MRD的輸出阻尼可以在外加磁場（一般通過電流大小來控制）作用下實現動態(tài)調節(jié)，并以其結構簡單、可調范圍寬、響應快、可靠性和適應性強等優(yōu)點被廣泛應用于建筑［2-4］、汽車［5-8］、航空［9-10］、航海［11-12］等領域。

一般阻尼器多應用于結構的減震、隔振及減擺等維持系統(tǒng)的穩(wěn)定性方面。在削弱振動過程中往往由于非線性因素［13］（本文MRD的非線性因素主要是MRF）的存在，使系統(tǒng)的運動具有非線性特性，非線性特性對多體系統(tǒng)動力學建模［14-15］的影響是不可避免的。黎崛珉等［16］研究了剛度非線性與阻尼非線性對隔振系統(tǒng)的影響。姚文莉等［17］將非理想系統(tǒng)的動力學建模問題描述為目標函數為高斯拘束函數、優(yōu)化變量為質點加速度的約束最優(yōu)化問題，并驗證其方法的有效性。Zhu［18，19］等人研究了帶有阻尼器的空氣彈簧非線性動力學模型，用來描述車輛減震系統(tǒng)以提高乘坐舒適性。

為了準確模擬MRD復雜的非線性阻尼力特性，實現對MRD精準控制，研究人員基于多體動力學［20］及MRF的材料特性對MRD的力學建模方法進行了深入探究［21-22］。Shou［23］、李趙春［24］等人研究了MRD在沖擊載荷作用下動力學響應特性模型。Ma［25］、Gurubasavaraju［26］、Yang［27］等研究了MRF剪切應力模型，對MRF電導率特性進行了 實 驗 研 究。Zolfagharian［28］、Peng［29］、陳 世嵬［30］等人對MRD進行參數化建模并對模型的參數辨識方法進行了研究。龔微［31］等提出了一種可以求解阻尼器逆模型的3段線性變阻尼恢復力模型，適用于MRD控制器的設計。梅真等［32］對MRD進行動力性能測試，建立了參數化模型和非參數化模型，通過對比2種建模方式的結果，發(fā)現非參數化模型的精度高于參數化模型。

以上研究結果表明，由于MRF的流變機理復雜，目前尚無統(tǒng)一、精確的力學模型描述，使得參數化模型的建模方式復雜化。非參數化模型是基于實驗測試結果建立的，利用數值計算方法對模型參數進行辨識，得到較為精確化的模型，可簡化建模方式和提高模型精度。孔祥東等［33］提出了Bingham-多項式力學模型，并與試驗測試結果進行對比，解決了Bingham模型精度低、多項式模型高階不穩(wěn)定的問題。祝世興等［34］提出用傅立葉級數來表達MRD的非線性關系，改進了線性模型精度差的問題。劉永強等［35］提出一種對MRD模型進行參數辨識的新方法，分析模型參數對MRD輸出阻尼的影響規(guī)律，獲取調控阻尼力的最佳參數，并擬合參數值和電流之間的函數關系。胡海剛等［36］提出一種將遺傳算法和最小二乘法相結合的方法，對MRD的力學模型進行參數辨識，并用Simulink工具箱進行仿真驗證，所得結果精度較高，可靠性滿足要求。周勇等［37］提出一種改進型RBF神經網絡模型，基于實驗實測結果，利用線性插值法建立模型算法，并與改進前的模型進行對比，發(fā)現改進后模型精度更高，可以更加可靠地描述MRD的動力學特性。

現有的非參數化模型參數較多，建模方式較為復雜，分段模型存在分段方式模糊、拐點不明確等問題。為了精確描述磁流變阻尼器在不同電流下力與速度之間的非線性關系、減少力學模型參數、簡化分段模型的分段方式，本文作者通過對臺架試驗的實驗結果進行分析，建立Rational-Bingham模型，采用最小二乘法對模型參數進行辨識，基于臺架試驗結果對Rational-Bingham模型和多項式-Bingham模型的精度進行對比研究。

1 實驗

實驗所采用的磁流變阻尼器結構原理圖如圖1，其中缸筒內填充MRF。

圖1 磁流變阻尼器結構原理圖

MRD實物圖如圖2，通過活塞桿引出的導線可以為MRD施加不同大小的電流。阻尼特性測試平臺如圖3所示，將MRD的上端固定，下端活塞端連接激振臺中心，激振臺通過液壓作動為阻尼器提供正弦位移激勵。

實驗采用正弦位移激勵，以振幅為2 mm、頻率為3 Hz為例，通過圖2引出的導線分別對MRD施加0、0.5、1.0、1.5、2.0、2.5 A的電流，研究在不同磁場強度下MRD的阻尼特性。將不同電流下力-位移關系轉換為力-速度關系，結果如圖4所示，其中橫坐標表示MRD的輸入速度v，縱坐標表示MRD的輸出阻尼力F。

圖2 實物圖

圖3 阻尼特性測試平臺

圖4 不同電流下MRD力-速度關系

從不同電流下MRD的動態(tài)響應曲線可以看出，將速度換向處（A處和C處）和加速度換向處（B處和D處）作為分段點進行分段建?？蓪⒁粋€運動周期分為2段線性曲線和2段為非線性曲線。

2 Bingham模型中MRF力學特性

MRF在無外加磁場的影響下可近似為牛頓流體，此時MRF的剪切應力與剪切應變率之間的關系為：

式中：τ為剪切應力；η為黏度系數；γ為剪切應變率。

MRF在外加磁場影響下具有黏塑性流體特性，如圖5所示，在不同的磁場強度下，流動模式可分為屈服前區(qū)、屈服區(qū)、屈服后區(qū)等3種變化過程。其中MRF在屈服前區(qū)剪切應力隨剪切應變率的變化呈正相關；在屈服區(qū)剪切應力隨剪切應變率的變化呈不穩(wěn)定狀態(tài)；在屈服后區(qū)呈穩(wěn)態(tài)，剪切應力不隨剪切應變率變化。

圖5 MRF的剪切應力與剪切應變率關系

出于本文的實際工況及MRD的工作狀態(tài)考慮，選取穩(wěn)態(tài)工作狀態(tài)下MRD中MRF的所處狀態(tài)進行研究，此狀態(tài)下MRF穩(wěn)定流動處于屈服后區(qū)，剪切應力為定值，MRF黏度為定值。

Bingham模型將MRF的剪切應力描述為：

式中：τy為MRF的剪切應力；τ0為MRF的臨界剪切應力；H為MRF所處磁場強度；γ為MRF的剪切應變率；η為MRF達到屈服后區(qū)的黏度系數。

Bingham模型中MRF在不同磁場強度影響下剪切應力與剪切應變率的量綱為一化關系如圖6所示。

圖6 Bingham模型中MRF本構關系

由圖6可知，磁場強度的變化會改變Bingham模型中的臨界剪切應力，而MRF在剪切屈服前后的黏度η在Bingham模型中被視為定值。

描述MRD輸出力與激勵速度之間線性關系的常用 Bingham力學模型，是由 Stanway等［38］提出的由庫倫摩擦力和黏滯阻尼器構成的理想化模型，即：

式中：fk為庫倫阻尼力，與MRF的屈服應力相關；x為MRD的輸出位移；˙x為MRD的輸出速度；c0為黏滯阻尼系數；f0為補償元件存在而產生的力，結合本文MRD實際結構及實際工況此項可忽略。

將圖4的實驗結果與Bingham力學模型進行對比分析，可以發(fā)現：圖4的線性部分在理論上近似表示MRF達到屈服后期且黏度不變時MRD的運動狀態(tài)，可通過Bingham模型模擬，但其非線性區(qū)域的滯回特性無法用Bingham模型進行精確描述。本文針對Bingham模型對描述MRD非線性區(qū)域力學特性存在的局限性問題，對傳統(tǒng)的Bingham模型進行改進，從而獲得更準確的磁流變阻尼器非線性力學模型。

3 力學建模與參數辨識

本文在Bingham模型基礎上進行分段建模，對圖4的線性部分采用Bingham模型，非線性部分分別采用多項式函數和Rational函數進行描述。

3.1 多項式-Bingham模型建立與參數辨識

對圖4進行分析可知，B-C段和D-A段所表示的MRD輸出力與速度之間的線性關系可通過Bingham模型進行描述，為了在形式上進一步簡化分段模型，本文將式（3）中fksgn（˙x）視為一個整體 f，c0視為 c，則參數 f和 c為待識別參數。A-B段和C-D段為非線性部分，首先采用多項式函數進行描述。為了防止高階振蕩失真和更好地描述滯回特性，擬采用3次多項式建模，可得多項式-Bingham模型（以下簡稱為多項式函數模型），其表達式為：

式中：v為MRD的速度；˙v為MRD的加速度；c1、f1、p11、p12…為待識別參數；式（4）中函數分別依次對應表示圖4中D-A段、A-B段、B-C段、CD段。

本文采用最小二程曲線擬合法對模型參數進行辨識，最小二成曲線擬合法的基本原理是對給定數據（xi，yi）（i＝0，1，…，n），在函數類 φ中，求得g（x）∈φ，使誤差的平方和R最小，可表示為：

從幾何意義上來說，就是找到與給定點（xi，yi）（i＝0，1，…，n）的距離平方和為最小的函數曲線 y＝g（x），函數 g（x）稱為最小二乘解或擬合函數，求擬合函數g（x）的方法稱為最小二乘曲線擬合法。

在振幅為2 mm、頻率為3 Hz的振動工況下，將加載電流分別為 0、0.5、1.0、1.5、2.0、2.5 A的實驗結果分別對多項式函數模型進行參數辨識，得到不同電流下對應的參數值，并建立不同參數與電流之間的關系。具體過程以加載電流為1.5 A的參數辨識結果為例，利用Matlab擬合工具箱對 CD段系數為 p21、p22、p23、p24的多項式函數進行曲線擬合，擬合結果如圖7所示，可得CD段的各參數 p21、p22、p23、p24分別為：0.000 016、0.001 25、0.035 35、0.063 13。

根據擬合結果可以看出，多項式函數與實際散點之間的擬合偏差較小，因此三階多項式具有一定的有效性。通過對多項式函數模型進行參數辨識，可得到不同電流下多項式函數模型的各參數值，部分結果如表1所示。

圖7 1.5 A電流下C-D段多項式函數的擬合結果

表1 多項式-Bingham模型參數辨識結果

通過對參數辨識結果進行分析可知，多項式函數模型各參數與電流之間為非線性關系，作者擬采用三階多項式來建立各個參數與電流之間的關系，可得到各參數與電流的函數關系，部分結果如下：

以上為多項式函數模型中各參數與電流之間建立的映射關系，結合式（4），可以得到完整的多項式函數模型，其中模型的輸入分別為電流和速度，輸出為對應輸入電流和速度下的輸出阻尼力。

3.2 Rational-Bingham模型建立與參數辨識

為了更好地描述同一振動工況不同電流下MRD的動態(tài)響應，本文建立了有理函數模型與多項式函數模型表進行對比研究。本質上有理函數（rational函數）是由2個多項式的商表示的函數，表示如下：

式中：P（x）為分子多項式函數；Q（x）為分母多項式函數。

當多項式擬合效果不夠好時，Rational函數擬合在精度上有一定的優(yōu)越性［39］。為了減少模型參數，將Rational函數的分子分母擬采用一次線性函數進行建模，可得到Rational-Bingham模型（以下簡稱為Rational函數模型）為：

式中：v為MRD的速度；˙v為MRD的加速度；a1、b1、d1、a2、b2、d2為待識別參數；式中函數分別對應表示圖4中 D-A段、A-B段、B-C段、CD段。

同樣在振幅為2 mm、頻率為3 Hz的振動工況下，將加載電流分別為 0、0.5、1.0、1.5、2.0、2.5 A實驗結果分別對Rational函數模型進行參數辨識，得到不同電流下對應的參數值，并建立不同參數與電流之間的關系。具體過程仍以加載電流為1.5 A的參數辨識結果為例，利用Matlab中的擬合工具箱對CD段系數為a2、b2、d2的Rational函數進行曲線擬合，擬合結果如圖8所示，可得CD段的各參數 a2、b2、d2分別為：-0.591 7、0.677 1、-11.78。

圖8 1.5 A電流下C-D段Rational函數的擬合結果

根據擬合結果可以看出，與多項式函數模型相比，Rational函數模型與實際散點的擬合效果好一些，通過對Rational函數模型進行參數辨識，得到不同電流下Rational函數模型中各參數值，部分結果如表2所示。

表2 Rational-Bingham模型參數辨識結果

同樣，對Rational函數模型參數辨識結果進行分析可知，各參數與電流為非線性關系，選擇采用三階多項式來建立各參數與電流之間的關系，得到各參數與電流之間的函數關系，部分結果如下：

以上為Rational函數模型中各參數與電流之間建立的映射關系，結合式（8）可得到完整的Rational函數模型，與多項式函數模型相似，Rational函數模型的輸入分別為電流和速度，輸出為對應輸入電流和速度下的輸出阻尼力。

4 模型仿真分析

本文在Matlab／Simulink環(huán)境中對模型在不同工況下力與速度之間的非線性關系進行仿真分析。如圖9所示，模型輸入分別為可以自定義的輸入電流及正弦位移激勵，模型輸出為MRD輸出阻尼力與速度之間的關系。其中SubSystem模塊中分別為多項式函數模型和Rational函數模型，是仿真系統(tǒng)的子系統(tǒng)。系統(tǒng)可分別輸出在不同工況下多項式函數模型和Rational函數模型的F-v仿真結果。

圖9 模型的Simulink仿真框圖

4.1 辨識工況（電流不同）

在振幅A為2 mm、頻率f為3 Hz的振動工況下，將不同電流分別輸入多項式函數模型仿真系統(tǒng)和Rational函數模型仿真系統(tǒng)，可分別得到輸出F-v曲線并與實驗結果生成對比，受篇幅限制，本節(jié)以施加電流大小分別為0、0.5、1.0、1.5 A為例，所得結果如圖10所示。

圖10 辨識工況下兩模型仿真結果與實驗結果

通過分析圖10中模型仿真結果與實驗結果對比，可以得出以下結果：

1）從整體上看，在同一振動工況不同電流下2種模型的擬合效果都比較好，都可以不同程度較為準確地描述經典Bingham力學模型無法模擬的非線性滯回區(qū)域。

2）從局部上看，多項式函數模型在低于0.5 A時高速段存在局部不穩(wěn)定性、在速度換向處存在擬合效果不佳的缺陷，而Rational函數模型可以修正多項式函數這一缺陷，在穩(wěn)定性上更加符合實際。

4.2 驗證工況（激勵頻率、振幅不同）

為了分別檢驗2種模型在同一電流不同振動工況下的適用性，將圖9中輸入電流設為I＝1.0 A固定不變，首先改變輸入正弦位移激勵信號的頻率，以振幅A為2 mm、頻率f為5 Hz為例，所得多項式函數模型和Rational函數模型的仿真結果與實驗結果對比如圖11所示。

圖11 I＝1.0 A，A＝2 mm，f＝5 Hz工況下，兩模型仿真結果與實驗結果

電流I＝1.0 A固定不變，改變輸入位移激勵信號的振幅，以振幅A為4 mm、頻率f為3 Hz的振動工況為例，所得多項式函數模型和Rational函數模型的仿真結果與實驗結果如圖12所示。

通過分析改變正弦激勵輸入信號的頻率或振幅后，2種模型的仿真結果與實驗結果對比（見圖11、圖12），可以得出以下結果：

1）多項式函數模型存在的缺陷更加顯著，改變頻率或振幅后，多項式函數模型在一個行程周期的高速段均存在不收斂的現象。線性區(qū)的仿真結果趨近于實驗結果，非線性低速區(qū)的仿真結果存在一定程度的誤差。

圖12 I＝1.0 A，A＝4 mm，f＝3 Hz工況下，兩模型仿真結果與實驗結果

2）為解決多項式函數模型在改變振動工況后高速段存在不收斂問題，可將多項式函數與線性函數相交之后的高速段取線性函數，低速非線性區(qū)用多項式函數來表達，這將改變分段模型的分段方式，因此多項式函數模型在本文的分段方式下對同一電流不同振動工況的適用性較弱。

3）在改變頻率或振幅后，Rational函數模型在非線性區(qū)的表達存在不同程度的誤差，但總體仿真結果與實驗結果相符，改變振動工況后仍然可以較好地模擬MRD輸出阻尼力的最大值。此外Rational函數模型可以修正多項式函數模型高速段不收斂的震蕩現象，對同一輸入電流不同振動工況下的適用性較強。

2種模型的非線性區(qū)均存在一定程度誤差的原因是：不同振動工況（正弦輸入信號）直接影響模型的輸入速度，所以簡化模型中各參數與正弦輸入信號的關系，在同一振動工況下建立各參數與電流之間的絕對關系，使得改變振動工況后非線性低速區(qū)的仿真結果存在一定程度的誤差。

5 模型誤差對比評估

2種非線性模型與實驗結果之間的誤差可通過均方根誤差RMSE來進行評價，計算公式如下：

式中：RMSE為均方根誤差；yi為擬合數據為實驗數據。

5.1 辨識工況

對同一振動工況（A＝2 mm，f＝3 Hz）不同電流下，2種模型的仿真結果與實驗結果各段的均方根誤差值及Rational函數與多項式函數相比Rational函數的誤差降低率，如表3所示。

CRMSE（累積均方根誤差）計算公式表示為：

ERR（誤差降低率）計算公式表示為：

表3 2種模型在不同電流下的均方根誤差

通過對表3進行分析，可以得出以下結果：

1）2種模型在同一振動工況不同電流下各段的擬合精度平均可達0.99以上；

2）不同電流下，Rational函數模型的累積擬合誤差與多項式函數模型的累積擬合誤差相比均有一定程度的降低，與2種模型各自的仿真結果一致；

3）在電流低于0.5 A時，Rational函數模型的誤差降低率隨著電流的增大而增大，電流高于0.5 A后Rational函數模型的誤差降低率隨著電流的增大而減小。

5.2 驗證工況

對同一電流（I＝1.0 A）不同振動工況（分別改變正弦位移激勵輸入信號的頻率和振幅）下，2種模型的仿真結果與實驗結果各段的均方根誤差值及Rational函數模型與多項式函數模型相比Rational函數模型的誤差降低率，如表4所示。

通過對表4進行分析，可以得出以下結果：

1）2種模型線性區(qū)的誤差結果較為接近，與模型結構相符合；非線性區(qū)多項式函數模型的誤差與Rational函數模型相比均有不同程度的增大，與仿真結果一致；

2）同一電流不同振動工況下Rational函數模型與多項式函數模型相比Rational函數模型的誤差降低率均在20%以上，結合圖11、圖12的仿真結果可知，Rational函數模型的適用性強于多項式函數模型。

表4 2種模型在不同振動工況下的均方根誤差

6 結論

1）相比于多項式函數模型，Rational函數模型的分子、分母項均為一次線性函數，模型參數較少，且模型結構簡單，分段方式明確。

2）在相同振動工況下，Rational函數模型可改進多項式函數模型在電流低于0.5 A時高速段存在的不穩(wěn)定性和在速度換向拐點處擬合效果不佳的現象，之后隨著電流增大，2種模型擬合誤差的差距逐漸縮小。Rational函數模型整體上可更加精確地描述MRD在不同電流下的非線性動態(tài)特性。

3）在相同電流下，改變輸入正弦位移激勵的頻率或振幅，通過不同頻率、不同振幅激勵工況的驗證表明，Rational函數模型與多項式函數模型相比，具有更高的精度和更好的適用性。