部分充液罐車內液體晃動的瞬態響應分析

王瓊瑤， 蔣開洪， Subhash Rakheja, 上官文斌

(1.華南理工大學 機械與汽車工程學院，廣州 510640；2.寧波拓普集團股份有限公司，寧波 315800)

隨著我國物流業和交通運輸業的迅猛發展，液罐車已成為石油、液化天然氣、危險化學品運輸的主要交通工具。由于液體貨物密度的不同以及最大載荷限制等原因，液罐車通常處于部分充液狀態。當罐車在制動或者轉向時，液體貨物未受限制的自由面將會發生顯著的晃動現象。液體的大幅晃動將會改變車軸的載荷分布，產生附加的力和力矩，從而對罐車的行駛性能造成嚴重的不利影響[1]。事故統計表明，相比較于其他的載貨車輛，液罐車更容易發生交通事故，且造成嚴重的人員傷亡和財產損失[2]。因此，研究液罐車內液體的晃動具有重要的現實意義。

對液體晃動現象的研究可追溯到1960年，然而研究主要集中在對航天器、船舶以及地面固定容器內液體晃動進行分析。對公路罐體內液體晃動的研究方法主要分為三類：準靜態方法，等效力學模型法以及液體晃動動力學模型法。準靜態方法主要用來預測液體晃動達到穩態時的動力學響應，且該方法只適用于不帶防波板的罐體[3]。然而，研究表明，液體晃動瞬態時的晃動力明顯大于穩態時的作用力[4]。等效力學模型可以用來分析液體的瞬態晃動，而且模型可以很方便地與車輛模型進行耦合[5]。包光偉[6]以油罐車為背景，采用Galerkin方法建立了貯箱平動激勵下離散的液體受迫晃動方程,建立了三維受迫晃動響應的等效力學模型。結果表明,縱向晃動固有頻率隨貯箱長度而下降,晃動質量隨貯箱長度而增加。等效力學模型的不足之處是模型參數需要從線性晃動理論或者試驗獲得，而且只限于對不帶防波板罐體內的液體晃動進行分析。因此，盡管防波板是一種有效的抑制液體晃動的裝置，但準靜態方法和等效力學模型法均不能分析防波板在抑制液體晃動方面的影響。液體晃動動力學模型法又可以細分為線性動力學方法以及CFD方法。這兩種方法均可對帶/不帶防波板的罐體內的液體晃動行為進行分析。線性動力學方法只限于分析液體的小幅度晃動，且需假設液體為不可壓縮，無黏性，無旋轉的流體。計算流體動力學(CFD)方法一般用于分析液體的大幅度晃動，由于該方法需要大量的計算量，關于其與車輛動力學模型耦合分析的研究較少[7-8]。

橫向防波板一般用來限制液罐車內液體貨物的晃動，尤其是在罐車制動的過程中。這類防波板通常中心開有一個圓孔，底部開有一個近似半圓形的通油孔。然而，用線性晃動理論和CFD方法對液體在帶防波板的罐體內的晃動進行分析的研究較少。Kolaei等[9]運用邊界元方法分析了水平圓柱形罐體內液體的晃動現象，研究并考慮了不同類型的縱向防波板。研究結果表明，固定于罐體頂部的部分防波板在常見的充液比的情況下能最有效地抑制液體的晃動，固定于罐體底部的防波板只在充液比較少的情況下才有效。劉小民等[10]以帶有單個防波板的運動罐體模型為基礎，對運動罐車在剎車過程中罐體內液體的晃動現象進行了數值分析，研究結果表明：罐體端面受力隨著充液比的增加而增大，當充液比(液體體積與罐體總體積的比值)為0.8時，罐體的總體受力最大。然而，以上研究中的防波板板均為平板，未考慮防波板的曲率對液體晃動的影響，也未考慮防波板的其他幾何參數的影響，比如防波板的傾斜角，開孔的形狀、位置和大小等。

本文建立了液罐車罐體內液體晃動的三維CFD模型，并應用此模型研究了防波板的幾何參數(曲率，開孔的大小、形狀及位置，防波板的傾斜角等)對液體瞬態晃動時的載荷轉移及晃動力的影響。模型的有效性通過與準靜態模型以及試驗結果的對比分析進行了驗證。用橫向和縱向載荷轉移量，作用在罐體端面、防波板上的力的穩態值和峰值以及晃動產生的俯仰力矩和側傾力矩評估了防波板的設計參數的影響。此外，對防波板與罐體端面以及相鄰防波板間的空氣壓力形成機理進行了分析，并研究了空氣壓力對載荷轉移以及晃動力產生的影響。

1 液罐車液體晃動的動力學模型

本文通過FLUENT求解器對液體晃動的動力學問題進行分析，如圖1所示。數值計算采用兩相流模型，即液相和氣相。氣相以空氣為研究對象，由于液體晃動過程中會出現氣體壓縮的現象，空氣采用可壓縮的理想氣體。氣液兩相密度差異較大，忽略了液體黏性力的影響。液體與罐體壁面、防波板間的相互作用可能會引起液體的大幅晃動甚至是液體飛濺現象，故流體域采用湍流計算模型。在氣液兩相流的數值計算中，兩相之間存在隨時間變化的運動界面，界面呈現多種多樣的形狀，本文采用的流體體積法(Volume of Fluid, VOF)能夠很好地分析液體自由面的變形，有效追蹤液體自由面的瞬態位置

▽·(uλ)=0

(1)

式中：u是流體速度矢量；▽是梯度算子；λ是流體體積分數。當λ等于0或者1時，表示網格單元分別為氣相和液相占據，當λ介于0和1之間時，則表示氣相和液相的交界面。

圖1 帶橫向防波板的液體晃動模型Fig.1 The liquid slosh model in a tank with baffles

給定流場Z方向重力加速度以模擬重力場。給定流場X方向和Y方向一定的加速度以分別模擬油罐車在制動和轉向時的工況，加速度隨時間的變化關系如下

(2)

式中：g(t)表示加速度；δ表示上升時間；G表示加速度的穩態值；β是使加速度值從上升期到平穩期平滑過渡的系數，取0.2；Ta是過渡時間。加速度曲線如圖2所示。

圖2 圓滑過渡的斜坡階躍加速度激勵Fig.2 Rounded ramp-step acceleration excitation

液體晃動時會使得液體的重心實時發生變化，液體的重心坐標(Xcg,Ycg,Zcg) 相對于罐體幾何中心的變化可由下式表示

(3)

記液體晃動時液體與罐體壁面接觸的區域為A1，空氣與罐體壁面接觸的區域為A2，如圖3所示。對于單元e，其對罐體壁面產生的作用力可由單元中心的壓強Pe與單元的面積矢量Ae求出。因而，液體晃動產生的合力為

圖3 液體的載荷轉移量，力以及力矩的計算原理圖Fig.3 Schematic diagram of computation of load shift, force and moment

(4)

式中：Fi為晃動合力沿i(i=x,y,z)軸的分量；Aei為單元e的面積矢量沿i軸的投影面積。俯仰力矩，側偏力矩以及側傾力矩的參考點均為罐體幾何中心點在罐體底部的投影點O′,通過單個單元的力矩在整個濕周的積分求得

(5)

式中：Fe為單元e產生的晃動力矢量；le為單元相對于點O′的位置矢量；M為力矩矢量。

液體晃動會對罐體壁面、端面和防波板產生作用力。因此，晃動產生的合力可以看成作用在罐體壁面的力Ft和作用在防波板上的力Fb的和，即

F=Ft+Fb

(6)

式(4)～(6)計算得到的力和力矩只考慮了液相部分，忽略了氣相(空氣)對罐體產生的作用力。當充液比較大時，防波板的孔將完全浸沒在液體中，如圖1(b)所示，在加速度激勵較大的情況下，防波板與罐體端面間或相鄰兩個防波板間空氣壓力將迅速增大。若考慮空氣壓力影響，液體晃動產生的總合力為

FT=F+Fa

(7)

式中：FT為總合力矢量；Fa為空氣產生的作用力矢量，其求解方式與液體晃動產生的作用力相似，為氣相單元e產生的作用力在區域A2的積分

(8)

在加速度激勵作用下，晃動力的穩態值與加速度的大小成線性關系。然而，晃動力的瞬態值遠大于其穩態值。所以在分析防波板在抑制液體晃動的作用時，有必要同時考慮力、力矩的瞬態值和穩態值。防波板的不同設計因素的作用可以通過瞬態放大因子來進行評估，其定義為作用在罐體和防波板上的力或力矩的最大值和相應的穩定值的比值，即

(9)

2 模型有效性檢驗

2.1 模型驗證方法

提出的CFD模型能夠有效計算液體在橫向或縱向加速度激勵作用下晃動產生的力和力矩。然而模型的有效性需要通過相關的理論或者試驗來驗證。常見的的驗證方法主要有準靜態模型法、頻率驗證法以及試驗驗證法。

2.1.1 準靜態模型法

在準靜態模型中，液體的自由面假設為平面，見圖4。

圖4 準靜態模型Fig.4 Quasi-static model

在給定的加速度激勵下，作用在罐體上的力和力矩為恒定值

(10)

2.1.2 頻率驗證法

液體自由晃動的頻率由罐體的形狀、大小以及充液比有關。Kobayashi等[11]1989年提出了圓柱形罐體在不同充液比時液體晃動頻率f的計算表達式

(11)

式中：g為重力加速度；L為液體自由面長度；h為靜態時液面的高度。在試驗中，給罐體一個橫向或者縱向簡諧加速度激勵，液罐內液體將會在加速度激勵的作用下晃動。當加速度激勵停止后，液體開始自由晃動；在CFD模型中，在流場的縱向和橫向分別給定一個持續時間為0.06 s的三角波脈沖加速度激勵，液體分別在俯仰平面和側傾平面自由晃動。通過對比試驗中和模型中液體自由晃動的頻率，即可進一步驗證模型的有效性。

2.1.3 試驗驗證法

為了進一步驗證模型的有效性，用0.264∶1的試驗模型(不帶防波板)進行驗證。試驗用的罐體由有機玻璃材料制成，罐體壁厚為10 mm，目的是盡量減小液體晃動過程中罐體壁面的振動。試驗是在振動試驗臺上進行。試驗的結構見圖5，罐體通過支撐架固定在上支撐板上，上支撐板通過四個傳感器與下支撐板連接，下支撐板固定在振動臺上。支撐板具體的結構尺寸見表1。

圖5 試驗模型結構圖Fig.5 The set-up diagram of the experimental model表1 支撐板的尺寸Tab.1 The dimension of the support plate

參數d/mw/mh/m數值1.580.70.03

試驗用的液體為水，充液比為50%，對罐體模型在縱向和橫向分別施加ax=0.5sin(2π×0.37t)以及ay=1.0sin(2π×0.95t)的簡諧激勵，罐體內的水在簡諧激勵的作用下來回晃動，見圖6。

(a) 瞬態:縱向晃動(b) 瞬態:橫向晃動

圖6 罐體模型試驗圖

Fig.6 The experiment diagram

將試驗測得的晃動力與力矩分別與仿真得到的晃動力與力矩進行對比，驗證模型有效性。

2.2 驗證結果及討論

圖7(a)為不帶防波板的罐體在gx=0.6g的斜坡階躍加速度激勵下以及充液比分別為50%、70%和90%時，由本文的CFD模型獲得的縱向力穩態值與準靜態模型獲得的縱向力穩態值的對比。

(a) CFD模型與準靜態模型獲得的縱向晃動力穩態值得對比

(b) CFD模型與式(12)獲得的自由晃動頻率的對比圖7 CFD模型驗證Fig.7 The validation of CFD model

由圖7(a)可知，兩者差別非常小，誤差主要由網格單元計算液體的總體積時出現的偏差造成的。圖7(b)為本文的CFD模型獲得的液體晃動自由頻率與式(12)獲得的結果的對比，由圖可知，當充液比為30%～90%時，縱向晃動的頻率范圍為0.16～0.26 Hz，橫向晃動的頻率范圍為0.56～0.81 Hz。縱向晃動頻率和橫向晃動頻率均隨著充液比的增大而增大。且CFD模型獲得的結果與理論值十分接近。

圖8(a)和(b)分別為縱向激勵和橫向激勵下，試驗結果與仿真結果的對比。由圖可知，在縱向激勵和橫向激勵的分別作用下，仿真值和試驗值都基本吻合。

在縱向晃動中，由仿真獲得的縱向力峰值比試驗值大，主要原因是在試驗中，液體自由面在晃動過程中出現分離(液體飛濺)現象，而在仿真中這一現象并未出現或液面分離程度較低。分離的液體下落到罐體內，見圖6(a)，使得試驗獲得的垂向力峰值出現雙峰現象。由于俯仰力矩不僅與縱向力和垂向力有關，還與液體的縱向載荷轉移量有關，雙峰現象在俯仰力矩中出現較少(僅在大約t=5 s時刻)；

在橫向晃動中，由試驗獲得的橫向力和垂向力峰值均比仿真獲得的峰值大，主要原因是在試驗中，罐體的壁面為有機玻璃面，在液體晃動過程中，罐體壁面也會發生相應的小幅振動，而在仿真中，罐體壁面假設為剛體，故不會出現罐體壁面振動的情況。由于側傾力矩不僅與橫向力和垂向力有關，還與液體的橫向載荷轉移量有關，試驗結果與仿真結果差異較小。

試驗驗證結果說明模型是有效的，可以有效預測三維液罐內液體晃動產生的力和力矩。

(a) 縱向激勵下的晃動響應對比

(b) 橫向激勵下的晃動響應對比圖8 試驗數據與仿真數據對比Fig.8 The comparison between the experimental data and simulation data

3 帶防波板的液罐模型的仿真分析

以實際的液罐車的罐體尺寸為依據，建立了1∶1相應的簡化幾何模型。簡化的幾何模型尺寸如圖9所示。罐體的橫截面為圓形，半徑為1 015 mm。罐體的圓柱部分沿X軸方向的長度為7 550 mm。罐體端面的曲率半徑為2 030 mm。三個防波板沿水平方向固定于罐體的中間位置。防波板的曲率半徑為2 030 mm, 孔的直徑為D, 底端通油孔的半徑r為144 mm，其面積占罐體橫截面積的1%。坐標原點O取在罐體的幾何中心。罐體內的液體為燃油，其密度為850 kg/m3。本文中除特殊說明外，防波板的中孔的形狀默認為圓形，中孔位于防波板的幾何中心，孔徑為884mm。仿真的初始條件為：罐體的充液比為50%或70%；給定流體域橫向(Y方向)或者縱向(X方向)斜坡-階躍加速度激勵(如圖2所示)，以及垂向(Z方向)的重力加速度激勵。邊界條件為：罐體的壁面以及防波板均為剛性壁面，且壁面為不可滑移邊界。

圖9 帶三防波板的液罐幾何圖Fig.9 Geometry of the tank with three baffles

3.1 網格和時間步長無關性檢驗

液體晃動的非線性相應受網格大小以及時間步大小的影響很大。選用三種不同大小的網格數來檢驗網格無關性，其對應的網格數分別為36 190,53 040以及207 148。選用了三種時間步長，分別為0.02 s，0.01 s和0.002 s。圖10為在加速度激勵為gx=0.3g，gy=0.25g，充液比為50%的情況下，縱向力和橫向力隨時間的變化曲線。由圖可知，當網格數由36 190增大到53 040時，縱向力和橫向力的響應曲線無論是峰值還是穩態值均會略微增大，而當網格數從53 040增大到207 148時，縱向力和橫向力的響應曲線無論是峰值還是穩態值均變化很小，對應的曲線幾乎重合，當網格數為53 040左右時，得到的結果可以接受。同理可得出時間步長為0.01 s較為合適，此處不再贅述。

3.2 罐體內空氣壓力的變化

在不帶防波板的罐體內，液體晃動過程中，空氣壓力產生的作用力通常非常小，可以忽略[12]。然而在帶防波板的罐體內，當充液比較大防波板的孔完全浸沒在液體中時，防波板與罐體端面間或相鄰兩個防波板之間的空氣壓力將會非常大，如圖1(b)所示。例如，在防波板1的左右兩側的臨近區域分別取一個空氣單元，如圖11(a)所示，當加速度激勵為gx=0.6g，gy=0.25g，充液比為70%時，監測兩個空氣單元在液體晃動過程中的氣壓隨時間的變化，見圖11(b)。由圖11可知，晃動過程中，防波板1左右兩側的氣壓均發生了大幅震蕩，且左側氣壓值始終大于右側氣壓值，表明左側空間空氣的壓縮程度相對更高。左側氣壓的峰值約為27 kPa，穩態值約為16 kPa，氣壓變化頻率為1.85 Hz，遠大于液體晃動頻率0.58 Hz。

(a) 橫向力

(b) 縱向力圖10 網格數對瞬態晃動力的影響Fig.10 Effect of grid size on transient slosh force

(a) 防波板1左右兩側的空氣單元

(b) 防波板1兩側空氣單元氣壓變化的時間歷程圖11 防波板1兩側的空氣單元及其氣壓變化的時間歷程Fig.11 The air element on both sides of baffle 1 and the time history of air pressure of the two elements

罐體端面與防波板之間以及相鄰防波板之間氣壓的大幅變化將會影響作用在單個防波板上的合力。圖12為當加速度激勵為gx=0.6g，gy=0.25g，充液比為70%時，作用在防波板1上的縱向力合力以及液體邊界單元、氣體邊界單元分別作用在防波板1上的縱向力。由圖12可知，由氣壓產生的作用在防波板1上的力甚至比液體晃動作用在防波板1上的力還要大，且兩者方向相反，從而減小作用在防波板上的合力，表明氣壓的影響可以限制液體的晃動。

圖12 分別考慮氣相和液相時防波板1上承受的縱向力Fig.12 The longitudinal force exerted on baffle 1 considering air phase and liquid phase, respectively

3.3 防波板的設計參數及其對液體晃動的影響分析

3.3.1 防波板的設計參數

研究了防波板的不同設計參數對液體晃動的影響，這些設計參數主要有：防波板中孔的直徑大小D(圖13(a))、防波板中孔的位置高度H(圖13(b))以及防波板中孔的形狀S(圖13(c))等，如表2所示。

表2 防波板的設計參數Tab.2 The design parameters of baffles

表2中的粗體部分表示防波板設計參數中的參考選項，為了研究這些設計參數對液體晃動的影響，固定其中兩項，變化第三項。與表2中設計參數對應的防波板的幾何形狀見圖13，其中ξ為防波板的曲率。

(a) 中孔的直徑大小D(b) 孔的位置高度H

(c) 孔的形狀

圖13 防波板的設計參數

Fig.13 The design parameter of baffles

3.3.2 防波板的設計參數對液體晃動的影響

研究表明，液體在縱向加速度激勵下，其載荷轉移量以及晃動力與防波板的開孔面積關系密切，但這些研究涉及的防波板均是平板。本文研究了彎曲型防波板的開孔面積的大小對液體晃動響應的影響。由于開孔面積對橫向載荷轉移以及橫向晃動力影響較小，討論僅限于縱向載荷轉移和縱向力。

圖14(a)為在加速度激勵為gx=0.6g，gy=0.25g，充液比為50%和70%時，縱向力放大因子隨防波板圓孔直徑變化的關系。防波板的圓孔直徑分別為D=642 mm，780 mm以及884 mm，對應的防波板的開孔面積(包括通油孔面積)分別占罐體橫截面積的10.5%，15.8%和20%。由圖可知，隨著防波板的開孔面積增大，縱向力放大因子逐漸增大，與充液比的大小無關。主要原因是：隨著防波板的開孔面積的增大，縱向載荷轉移量增大，如圖14(b)所示，表明液體晃動的平均速度增大，使得晃動的劇烈程度增大，從而增大了縱向力的大小。

(a) 縱向力放大因子

(b) 縱向力載荷轉移量圖14 防波板的開孔面積對縱向晃動響應的影響Fig.14 The effect of opening area of baffles on liquid slosh response

除了防波板的開孔面積，防波板孔的形狀也會對液體晃動產生影響。圖15為在加速度激勵為gx=0.6g，gy=0.25g，充液比為50%和70%時，防波板的開孔形狀對縱向力放大因子以縱向力載荷轉移量的影響。孔的三種形狀分別為圓形，橢圓形以及十字架形，它們的開孔面積相等，與底部的通油孔一起均占罐體橫截面積的20%。由圖15(a)可知，與圓形孔和橢圓形孔的防波板相比，防波板孔為十字架形時，縱向力的峰值最小，與充液比無關。主要原因是：在十字架形孔周圍會形成渦流，渦流的形成會減小液體晃動的動能[13]，從而減小晃動的劇烈程度。然而，當充液比為70%時，十字架形孔的防波板將使液體縱向載荷轉移量更大，見圖15(b)。主要原因是：十字架形孔的上沿位置相比于圓形孔和橢圓形孔的上沿位置更高，如圖15(b)中的小圖所示。在中等或高充液比的情況下，比較不容易形成封閉空間，空氣壓力的變化相對較小，其對液體晃動的抑制作用相對較小。橢圓形孔由于其孔的上沿位置最低，容易形成封閉空間，液體晃動過程中，空氣壓力的變化較大，其對晃動的抑制作用也因此變大。因此對橢圓形孔的防波板而言，其引起的縱向載荷轉移量最小，如圖15(b)所示。

(a) 縱向力放大因子

(b) 縱向力載荷轉移量圖15 防波板的開孔形狀對縱向晃動響應的影響Fig.15 The effect of opening shape of baffles on liquid slosh response

受防波板孔的形狀對液體晃動影響的啟發，研究了防波板孔的位置高度對液體晃動的影響。圖16為在加速度激勵為gx=0.6g，gy=0.25g，充液比為50%和70%時，防波板孔位置高度對縱向力放大因子以及縱向載荷轉移的影響。由圖可知，當充液比為50%時，隨著孔的位置高度的增大，縱向力峰值減小，而縱向載荷轉移量卻隨之增大。當充液比為70%時，隨著孔的位置高度的增大，縱向力峰值以及縱向載荷轉移量均隨之增大。主要是因為孔的位置高度將影響防波板的有效面積以及罐體內空氣壓力的形成。當充液比為50%時，孔的位置高度增大，防波板的有效面積增大，因此縱向力峰值減小，同時，孔的位置較高時，罐體內防波板與罐體端面或相鄰防波板間的封閉空間不易形成，空氣壓力的變化很小，縱向載荷轉移此時較大；當充液為70%時，若孔的位置高度為H=715 mm，則孔完全浸沒于液體中，罐體端面與相鄰的防波板之間的封閉空間在晃動的一開始便形成了，液體在加速度激勵的作用下向罐體端面聚集并壓縮封閉空間，封閉空間的氣壓會突然增大(如圖11所示)。增大的氣壓會阻礙液體進一步向端面聚集，從而減小液體的縱向載荷轉移量以及液體晃動的幅值。隨著孔的位置高度逐漸增大，孔逐漸露出液面，封閉空間在晃動初始階段并未形成，氣壓幾乎不變(一個標準大氣壓)且其對晃動無影響，直至液體向罐體端面聚集過程中封閉空間形成。因此，當充液比為70%時，隨著孔的高度的增大，縱向力峰值以及縱向載荷轉移量均隨之增大。

(a) 縱向力放大因子

(b) 縱向力載荷轉移量圖16 防波板孔的位置高度對縱向晃動響應的影響Fig.16 The effect of opening height of baffles on liquid slosh response

4 結 論

(1) 建立了部分充液罐車流體動力學模型，并通過準靜態模型和試驗方法等驗證了模型的有效性。

(2) 研究了防波板的幾何參數對液體晃動的影響。研究表明，防波板的孔的大小、形狀以及位置高度對液體晃動的影響很大。通過減小防波板孔的面積以及降低孔的位置高度，均可提高防波板抑制液體晃動的能力。

(3) 提出了當防波板的孔的位置高度較低或充液比較大時，防波板與罐體端面或者相鄰防波板間空氣壓力的形成機制，研究表明空氣壓力可以較小液體的載荷轉移從而降低液體晃動的劇烈程度。因此空氣壓力的形成能有效提高防波板抑制液體晃動的能力。