陳 雨 路力權 張軍昌 陳 軍 胡宸瑋 曹佳宇

(西北農(nóng)林科技大學機械與電子工程學院， 陜西楊凌 712100)

0 引言

高地隙自走式噴霧機具有高效、環(huán)保等優(yōu)勢，得到了廣泛應用[1-3]。但噴霧機運行工況復雜，地隙高，輪距相對于車身高度偏窄[4-5]，作業(yè)時藥液易發(fā)生晃動。液體易于流動的特性和噴霧機運行狀態(tài)的改變使得藥液在外力作用下沖擊藥罐壁面，產(chǎn)生變化的動態(tài)負載，致使作用于整機的力矩增加，極大地降低了整機穩(wěn)定性與作業(yè)質(zhì)量。因此，研究藥罐內(nèi)液體晃動的動力學特性，對于抑制整機側傾與俯仰運動，提高噴霧機穩(wěn)定性與噴藥效果具有實際意義[6]。

目前，車輛罐體內(nèi)液體晃動動力學的研究方法主要有流體動力學法[7-8]、數(shù)值解析法[9]、等效力學模型法[10]等。其中，等效力學模型法是指通過力學模型近似描述液體晃動，將流體力學問題轉化為機械運動問題。由于其更容易納入到系統(tǒng)穩(wěn)定性分析、控制器設計以及系統(tǒng)模擬中，因此得到了廣泛的研究。例如KOLAEI等[11]將非滿載罐體液體晃動等效力學模型與半掛液罐車五自由度動力學模型進行耦合分析，確定了液罐車側傾穩(wěn)定特性。NOKHBATOLFOGHAHAI等[12]、李金輝等[13]將液體晃動彈簧-質(zhì)量模型整車多體動力學模型結合，分析了液體晃動對整車穩(wěn)定性的影響。NICOLSEN等[14]和GROSSI等[15]通過集成基于拉格朗日連續(xù)介質(zhì)的液體晃動模型和卡車多體系統(tǒng)模型，研究了液體晃動對罐車動力學特性的影響。鄭雪蓮等[16-17]、孫文財?shù)萚18]和趙偉強等[19]將液罐內(nèi)液體晃動等效為單擺力學模型，分析不同工況條件下變充液比罐內(nèi)液體晃動所帶來的沖擊效應對于罐式半掛車整車側傾穩(wěn)定性的影響，并以此提出了基于差動制動的液罐車防側翻控制策略。

綜上所述，液體晃動等效力學模型在罐式道路運輸車輛的穩(wěn)定性分析與控制研究中已廣泛應用。與罐體車輛相比，噴霧機質(zhì)心高，藥罐以圓柱和橢圓柱為主，藥液與底盤質(zhì)量比大且作業(yè)時藥液質(zhì)量時變，對于果園噴霧機[20-22]及中小型植保機[23-25]，受限于罐體體積，藥罐內(nèi)未放置防波板，藥液晃動對其整車穩(wěn)定性和安全性影響更大。研究構建可準確描述噴霧機液體縱向晃動過程的等效力學模型，可探明整機平順性及穩(wěn)定性等影響因素，進而為開展噴霧機底盤結構改進、懸架設計、罐體形狀和防晃裝置優(yōu)化、整機平順穩(wěn)定控制等提供依據(jù)?；诖?，本文以噴霧機藥液縱向晃動為例，基于液體晃動動力學特性及機械模型等效準則，研究建立藥罐液體縱向晃動等效力學模型，并通過仿真與實驗方法驗證模型的準確性。在此基礎上，建立1/2噴霧機四自由度垂向動力學模型，分析討論液體晃動與藥液時變對噴霧機行駛平順性能的影響。

1 液體縱向晃動等效力學模型建立

1.1 液體縱向晃動特性分析

高地隙自走式噴霧機藥罐形狀多為圓柱體或橢圓柱體，其縱向截面為矩形。建立罐內(nèi)液體晃動坐標系如圖1所示。取自由液面中心為坐標原點o，x為噴霧機前進方向，z為垂直方向。圖中，2a為罐體長度，H為自由液面高度，Ω為罐內(nèi)流體域。

圖1 罐體縱向晃動坐標系Fig.1 Coordinate system of tank longitudinal shaking

當外界激勵較小時，罐內(nèi)液體處于小幅晃動。在如圖1所示的平面坐標系oxz中，設h(x,z,t)為自由液面波高函數(shù)，引入流體相對速度勢函數(shù)Φ(x,z,t)，由于小幅晃動下流體可看作是無旋的，則有

(1)

液體粒子在罐體內(nèi)表面處的法向速度滿足

(2)

(3)

同時，自由液面上的液體粒子在z方向上的運動應滿足

(4)

(5)

式中ρ——流體密度，kg/m3

p——液體壓力，N

g——重力加速度，m/s2

去掉式(5)非線性項，可得

(6)

則罐內(nèi)液體產(chǎn)生的總壓力可寫為

(7)

在自由表面上，液體滿足的動力學條件為

(8)

式(1)～(4)、(8)構成了罐內(nèi)液體受迫晃動的完整動力學方程。

在罐內(nèi)液體小幅晃動時，其一階模態(tài)起主導作用。基于此，在對液體晃動的控制方程求解時不考慮其多階模態(tài)的疊加，僅考慮一階模態(tài)，此時可將液體的速度勢函數(shù)Φ(x,z,t)展開為

(9)

顯然，式(9)滿足式(1)～(3)，將其代入式(4)、(8)中得

(10)

(11)

(12)

再將式(12)回代至式(11)中得

(13)

式中ω1——罐內(nèi)液體一階晃動自然頻率

將式(13)右端的-x展開為三角級數(shù)，得到

(14)

由Duhamel積分得到t=0條件下的方程(14)的解為

(15)

由式(15)可得液體的液動壓力為

(16)

(17)

(18)

1.2 等效力學模型構建與參數(shù)求解

目前，使用較多的等效模型為單擺模型和彈簧-質(zhì)量-阻尼模型[27]。這兩種模型都可以較為準確地描述罐內(nèi)液體往復沖擊的非線性運動，而后者在低頻小幅的振動下對液體運動的描述更為精確。因此，基于液體縱向晃動特性分析，建立如圖2所示等效力學模型。圖中，G為液體質(zhì)心，hg為質(zhì)心高度；m0為固定質(zhì)量塊質(zhì)量，h0為固定質(zhì)量塊高度；mk、kk、ck分別為等效質(zhì)量塊質(zhì)量、等效彈簧剛度以及等效阻尼系數(shù)，hk為等效質(zhì)量塊高度，φ為車身俯仰角。

圖2 縱向彈簧-質(zhì)量-阻尼等效模型Fig.2 Longitudinal spring-mass-damp equivalent model

(19)

進而，根據(jù)系統(tǒng)動力學理論，該系統(tǒng)對罐壁在水平方向上產(chǎn)生的力與力矩為[20]

(20)

(21)

式中ωk——k階液體晃動固有頻率

根據(jù)前文對液體晃動過程的分析，在求解模型參數(shù)時同樣只分析其一階振動模態(tài)。根據(jù)力學相似原則，罐內(nèi)實際液體沖擊產(chǎn)生的力與力矩與等效模型產(chǎn)生的力與力矩相等，進而根據(jù)方程兩邊時間函數(shù)的系數(shù)相等可求得等效模型各參數(shù)，結果為

(22)

式中M——罐內(nèi)液體總質(zhì)量，kg

m1——等效移動質(zhì)量塊質(zhì)量，kg

(23)

式中h1——等效移動質(zhì)量塊高度，m

(24)

(25)

(26)

本研究中噴霧機藥罐截面半徑600 mm，長1 000 mm。 為分析不同充液比下液體晃動沖擊作用，根據(jù)其幾何尺寸及式(22)～(26)，以0.1為變化步長，分別求得充液比0.1～0.9時縱向等效力學模型各參數(shù)如表1、2所示。

表1 縱向等效模型等效質(zhì)量和高度Tab.1 Equivalent mass and height of longitudinal equivalent model

表2 縱向等效模型等效剛度和等效阻尼系數(shù)Tab.2 Equivalent damping and equivalent stiffness of longitudinal equivalent model

1.3 等效力學數(shù)值解析模型建立

縱向等效模型在求解過程中采用數(shù)值解析法，通過力學等價原則求解得到模型參數(shù)。為后續(xù)進一步驗證模型的準確性，需結合求得的參數(shù)及液體晃動力矩表達式建立等效力學數(shù)值解析模型。

根據(jù)式(19)～(21)求得在藥罐沿縱向晃動時藥液產(chǎn)生的作用力矩為

(27)

n——液體晃動模態(tài)階數(shù)，取1

結合式(27)在Matlab/Simulink平臺中搭建數(shù)值仿真模型，如圖3所示。

圖3 Simulink等效力學模型系統(tǒng)Fig.3 Simulink equivalent mechanics model system

2 等效模型有效性驗證

將等效數(shù)值解析模型在Matlab/Simulink中的仿真結果與液體晃動在Fluent中的仿真結果進行對比，以驗證所建等效模型的準確性。

2.1 Fluent仿真模型建立與驗證

在Fluent流體模擬軟件中對罐內(nèi)藥液晃動產(chǎn)生的沖擊作用進行模擬。按照噴霧機藥罐實際尺寸建立藥罐模型并標定罐內(nèi)液體區(qū)域如圖4所示。設置充液比變化范圍為0.1～0.9，變化步長為0.1。

圖4 Fluent仿真液罐模型Fig.4 Fluent simulation liquid tank model

根據(jù)液體受迫晃動的作用特點設置仿真求解器為基于壓力的瞬態(tài)流動求解器，罐內(nèi)同時存在水與空氣的情況下，設置其物理模型為VOF多相流模型，湍流模型為標準k-ε模型[22-23]。通過用戶自定義函數(shù)UDF(Userdefined function)實現(xiàn)對藥罐初始激勵的加載[24]。仿真時間步長0.01 s，仿真步數(shù)500步。

設計并制作了液體晃動實驗臺架以驗證Fluent仿真的準確性，如圖5所示，在該平臺上可對按照噴霧機藥罐等比例制作的小型液罐進行一定頻率和幅值的簡諧激勵。

圖5 模擬液體晃動實驗臺架Fig.5 Experimental bench for simulating liquid shaking1.液罐 2.滑塊導軌 3.液罐承托板及固定托 4.高速攝像機 5.控制器顯示單元 6.調(diào)速電機 7.攝影燈 8.曲柄偏心輪

由于液罐模型為圓柱形，難以直接測得液晃實驗中罐內(nèi)液體對罐壁產(chǎn)生的沖擊力，因此，本實驗通過在液體晃動過程中監(jiān)測液面狀態(tài)的變化，將其與Fluent仿真中的實時液面狀態(tài)進行對比，以驗證其仿真結果的有效性。液面監(jiān)測采用OLYMPUS公司i-SPEED TR高速攝像機，分辨率為1 280像素×1 024像素，最高幀率達10 000 f/s，適用于絕大多數(shù)測試，同時該攝像機配合獨特的控制器顯示元件(Control and display unit，CDU)來使用。

與仿真實驗中初始條件給予的速度激勵不同，在液晃模擬實驗中給予罐體的是位移激勵

Sx=0.02πsin(8πt)

(28)

根據(jù)式(28)，調(diào)節(jié)電機轉速及偏心輪偏心距。將獲得的高速攝像結果與Fluent實時液面狀態(tài)對比，圖6為充液比α=0.5條件下，不同時刻對比結果，左側為實驗結果，右側為仿真結果。

圖6 不同時刻實驗結果與仿真結果對比Fig.6 Comparison of experimental results and simulation results

由圖6可知，在相同時刻下，F(xiàn)luent仿真結果與模擬液晃實驗結果高度相似。為量化兩者液面匹配程度，借用Matlab中的grabit插件分別提取一個運動周期內(nèi)不同時刻的液面曲線，如圖7所示。

圖7 實驗與仿真結果液面曲線Fig.7 Comparison of liquid level curves between experimental and simulation results

從圖7中可以看出，在同一時刻下兩者的液面曲線在變化趨勢上基本一致，但在液面上水波高度的峰值存在略微差異。而從直觀曲線中難以評價兩種結果的液面曲線是否具有較高的擬合度。因此，在評價Fluent結果與高速攝像結果曲線的擬合度時，用correl函數(shù)來表示，即

(29)

根據(jù)式(29)計算兩種結果下液面曲線的相關系數(shù)，結果如表3所示。

表3 不同充液比下各時刻的相關系數(shù)Tab.3 Correlation coefficient results at each moment under different liquid filling conditions

結合correl函數(shù)定義，兩組數(shù)據(jù)相關系數(shù)越接近1時，其相關性越強，且呈現(xiàn)正相關。從表3可知，F(xiàn)luent仿真液面結果與高速攝像結果兩種曲線的相關系數(shù)在不同充液比時，不同時刻下基本處于0.85以上，說明Fluent仿真能夠在一定誤差范圍內(nèi)反映液體晃動的真實狀態(tài)，可進一步提高等效模型驗證的合理性。同時相關系數(shù)計算結果中亦存在小于0.6，甚至小于0.2的情形，由于進行液晃實驗時，液體晃動產(chǎn)生的慣性力帶給電機一定的運動阻力，因此，實際的晃動激勵可能與期望的原始激勵存在一定偏差，造成實驗的運動周期相對Fluent仿真存在一定滯后。此外，由于高速攝像過程中液罐處于運動狀態(tài)，因此存在拍攝角度差異，同時在進行液面曲線提取過程中亦存在一定誤差。

2.2 等效模型準確性驗證

根據(jù)圖3數(shù)值解析仿真模型在Matlab/Simulink中進行仿真，其中各模塊初始值設置為表1、2中的等效模型參數(shù)。選擇Fixed-step固定步長模式的ode5求解器，設置仿真步長為0.001 s，仿真步數(shù)為5 000步。通過step模塊與sinewave模塊給予系統(tǒng)一定初始激勵，在此激勵的作用下將力學系統(tǒng)產(chǎn)生的作用力矩結果記錄并導出。同時，使用Fluent流體模擬軟件對罐內(nèi)藥液晃動產(chǎn)生的力矩進行仿真模擬。

為提高仿真實驗驗證的精確性，結合噴霧機實際作業(yè)工況，選擇將階躍和簡諧兩種初始激勵條件下的仿真結果進行對比。其中，階躍激勵為

(30)

簡諧激勵為

vx=0.16πcos(8πt)

(31)

圖8為階躍激勵下仿真結果對比，圖9為簡諧激勵下仿真結果對比。

圖8 階躍激勵下仿真結果對比Fig.8 Comparison of results under step excitation

圖9 簡諧激勵下仿真結果對比Fig.9 Comparison of results under simple excitation

由圖8、9可知，在不同初始激勵條件下Matlab/Simulink仿真中等效模型數(shù)值分析與Fluent流體仿真的藥液晃動作用力矩隨時間變化規(guī)律基本一致，說明所建立的縱向等效力學模型的準確性與合理性。通過分析圖中的曲線誤差可以得出，在藥罐近乎空載即充液比為0.1時，兩者誤差略大。由于藥罐接近空載時，罐內(nèi)液體晃動的非線性特性較強，兩者仿真結果存在一定的誤差。此外，在階躍激勵條件下，最大誤差出現(xiàn)在激勵結束的時刻，由于仿真實驗中階躍激勵施加時間短暫，此時系統(tǒng)處于非穩(wěn)定狀態(tài)，出現(xiàn)響應滯后的情形，增大了誤差。

3 等效模型應用

3.1 1/2噴霧機垂向動力學模型建立

基于上述建立的等效力學模型與車輛動力學模型耦合，建立綜合藥罐液體晃動和不考慮藥罐液體晃動的1/2噴霧機垂向四自由度動力學模型，通過仿真對比，探明罐內(nèi)藥液晃動對噴霧機行駛過程中整機行駛平順性的影響。

以“地面-輪胎-懸架-車身”構成的系統(tǒng)為對象，建立不考慮藥罐液體晃動的1/2噴霧機四自由度垂向動力學模型，如圖10所示。傳統(tǒng)模型建立過程中，假設噴霧機車身與藥罐為剛體，不考慮罐體內(nèi)液體晃動和車輛姿態(tài)引起的輪胎載荷轉移的影響。圖中zb為質(zhì)心垂向位移；zF、zR為前后簧載質(zhì)量的垂向位移；ztF、ztR為前后輪胎垂向位移；zwF、zwR為路面位移激勵；msc為半車簧載質(zhì)量，此時的簧載質(zhì)量可看作是將罐內(nèi)液體及車架視為同一剛體的簡化結果；Isc為車身轉動慣量；ksF、csF、ksR、csR分別為前后懸架等效剛度與等效阻尼系數(shù)；ktF、ksR為車輪等效剛度；lf、lr為車身質(zhì)心至前后軸在地面的投影距離。

圖10 不含液體晃動的噴霧機半車動力學模型Fig.10 Dynamics model of sprayer half-car without liquid sloshing

結合圖10，根據(jù)車身質(zhì)心處的垂向位移zb和俯仰角φ，動力學方程可表示為

(32)

(33)

(34)

(35)

當俯仰角φ較小時，近似有

zF=zb+lfφ

(36)

zR=zb-lrφ

(37)

結合式(27)、(32)～(37)，建立綜合液體縱向晃動的1/2噴霧機非線性垂向動力學模型，如圖11所示。

圖11 含液體晃動的噴霧機半車動力學模型Fig.11 Dynamic model of sprayer half-car with liquid sloshing

圖11中msc1為罐內(nèi)液體質(zhì)量，msc2為車架質(zhì)量，且有msc1+msc2=msc。此時，由于要分析液體晃動的影響，故將液體與車架質(zhì)量分開考慮。動力學方程表示為

(38)

(39)

(40)

3.2 Simulink動力學仿真

為探究罐內(nèi)液體晃動對噴霧機行駛穩(wěn)定性的影響，結合式(32)～(40)在Matlab/Simulink中搭建動力學仿真模型，如圖12所示。

圖12 Simulink半車動力學仿真模型Fig.12 Simulink half-car dynamics simulation model

仿真參數(shù)根據(jù)本團隊已設計的小型噴霧機設置，半車車架質(zhì)量msc2為10.3 kg，前后非簧載質(zhì)量mtF、mtR為15.8 kg，懸架等效剛度ksF、ksR為16 000 N/m，懸架等效阻尼系數(shù)csF、csR為500 N·s/m，車輪等效剛度ktF、ktR為20 000 N/m，前車軸與質(zhì)心投影距離lf為0.75 m，后車軸與質(zhì)心投影距離lr為0.5 m，噴霧機行駛速度v為1.4 m/s，其中，仿真過程中隨著罐內(nèi)充液比的不同，半車簧載質(zhì)量msc1、液體質(zhì)量msc2、車身轉動慣量Isc如表4所示。

表4 噴霧機仿真參數(shù)Tab.4 Simulation parameters of sprayer

仿真中以生成的隨機路面為初始激勵信號。結合噴霧機實際作業(yè)路況，選擇仿真的路面激勵模型為E級路面[28-29]。采用一階濾波帶白噪聲法建立路面激勵時域模型

(41)

式中w(t)——路面位移輸入，m

q(t)——限帶白噪聲

v——噴霧機行駛車速，m/s

nd0——路面空間截止頻率,Hz

考慮前后車輪的路面激勵存在一定的時間延遲，結合前后軸距及噴霧機行駛速度，生成車輪隨機路面激勵信號如圖13所示。

圖13 前后車輪隨機路面激勵Fig.13 Random road input for front and rear wheels

將噴霧機車身加速度與俯仰角速度作為信號輸出，并以此作為評價指標，對比有無液體晃動時的仿真結果，如圖14所示。

圖14 有無液體晃動時仿真結果對比Fig.14 Comparison of simulation results

由圖14可知，當噴霧機在隨機農(nóng)田路面行駛時，由于藥罐內(nèi)液體晃動力作用，其車身加速度與俯仰角速度增大，降低了噴霧機行駛的平順性。說明噴霧機行駛作業(yè)時，若將藥罐藥液視為剛體，忽略液體晃動的沖擊作用時，難以準確分析噴霧機整機運動狀態(tài)及車身姿態(tài)的變化。因此，本文所建立的液體縱向晃動等效力學模型可準確分析噴霧機整機行駛平順性、穩(wěn)定性。

為進一步分析噴霧機藥液質(zhì)量時變下液體晃動對整機性能影響，本文對不同充液比條件下的噴霧機垂向動力學特性進行對比分析。圖15所示為充液比0.1、0.5、0.9下，噴霧機車身加速度與俯仰角速度的變化情況，并將不同充液比下兩評價指標的最大值與均方根進行對比分析，如表5所示。

圖15 不同充液比下仿真結果對比Fig.15 Comparison of simulation results at different filling ratios

由圖15可知，隨著藥液質(zhì)量的增大，液體晃動對噴霧機行駛作業(yè)的影響逐漸增大，且整車在路面激勵下響應的滯后效應越明顯。由表5可知，在相同質(zhì)量下，若忽略液體晃動影響，將其視為簧載車身質(zhì)量的一部分時，整車的俯仰角與車身加速度顯著降低。對比是否考慮液晃影響下的數(shù)據(jù)結果，充液質(zhì)量越大，兩者的結果差值越大，進一步說明了液體晃動對整車行駛性能的影響。此外，當充液比為0.8時，整車的俯仰角與車身加速度達到最大值，此時的噴霧機行駛平順性較差。由于空載與滿載(即充液比為0和1)時，不存在罐內(nèi)液體晃動作用，因此，表5中上述兩種特殊情況下是否考慮液體晃動的仿真結果相同。通過以上分析，本文所建立的藥罐液體晃動等效模型，在揭示液體晃動作用對噴霧機行駛穩(wěn)定性影響的同時，也能反映藥液時變因素對其所造成的影響，為進一步研究噴霧機不同作業(yè)工況和藥液質(zhì)量時變下底盤結構優(yōu)化及平順穩(wěn)定分析提供實施途徑和參考方案。

表5 不同充液比仿真結果對比Tab.5 Comparison of simulation results of each filling ratio

4 結論

(1)建立了噴霧機藥罐內(nèi)液體縱向晃動的彈簧-質(zhì)量-阻尼等效力學模型，并對模型參數(shù)進行了求解。應用Fluent建立了罐內(nèi)液體晃動仿真模型，設計并進行了模擬液體晃動實驗，實驗結果表明，實際晃動液面曲線與Fluent仿真液面曲線相關系數(shù)大于0.85，說明了基于Fluent的液體晃動仿真模型的準確性。

(2)應用Matlab/Simulink軟件，建立等效力學數(shù)值解析模型，結合Fluent流體仿真結果，對比分析了充液比為0.1、0.5、0.9時液體縱向晃動作用在容器壁的力矩變化規(guī)律。對比結果表明，所建立的等效力學模型力矩變化與Fluent中的力矩變化過程一致，驗證了所建立的液體縱向晃動等效力學模型的正確性。

(3)建立了綜合藥罐液體晃動和不考慮藥罐液體晃動的半車噴霧機四自由度垂向動力學模型，通過仿真對比分析表明，藥罐內(nèi)液體晃動使車身加速度與俯仰角速度增大，降低了噴霧機行駛的平順性。通過分析不同充液比條件下的噴霧機垂向動力學特性表明，當充液比為0.8時，整車行駛平順性較差，在進行轉場或噴藥作業(yè)時應盡量避免充液比為0.8的情況。