液罐車內(nèi)彈性擋板對液體晃動(dòng)的響應(yīng)分析*

徐 凱，王瓊瑤，林國勉

（五邑大學(xué)智能制造學(xué)部，廣東 江門 529000）

0 引言

液體晃動(dòng)現(xiàn)象常出現(xiàn)在部分充液容器中，并且涉及到工程領(lǐng)域的各個(gè)方面，包括航空航天領(lǐng)域、陸運(yùn)領(lǐng)域、海洋工程領(lǐng)域以及儲(chǔ)罐等。液體的劇烈晃動(dòng)可能導(dǎo)致充液容器的結(jié)構(gòu)破壞、液罐車側(cè)翻事故的頻發(fā)以及航天器發(fā)射的失敗等。因此，研究如何抑制充液容器內(nèi)液體的晃動(dòng)對工程應(yīng)用具有十分重要的意義。

自20 世紀(jì)60年代開始，Abramson 等[1]比較全面地開展了常重力環(huán)境下各種形狀的擋板的液體晃動(dòng)實(shí)驗(yàn)研究。Hasheminejad 和Mohammadi[2]使用線性勢理論的二維動(dòng)力學(xué)方法研究了液體在無黏性和不可壓縮的條件下，水平圓柱擋板容器的3 種擋板安置方式對液體晃動(dòng)的影響。結(jié)果表明：容器在高充液比情況時(shí)，底部安裝垂直擋板不是抑制部分填充容器內(nèi)液體晃動(dòng)的最有效的手段。Yu和Chu[3]研究了部分充滿液罐車的液體晃動(dòng)的非線性運(yùn)動(dòng)，研究結(jié)果表明：在充液比在90%時(shí)，瞬時(shí)側(cè)向力最小。Yan 和Rakheja[4]研究了部分充液的液罐車在有、無擋板的情況下其直線制動(dòng)性能特征。研究結(jié)果顯示當(dāng)罐體內(nèi)有擋板時(shí)，液罐車的制動(dòng)性能明顯提高。Wang 等[5]基于VOF （Volume of Fluid） 模型研究了無擋板，傳統(tǒng)剛性擋板和組合擋板在橫、縱方向上對液罐車運(yùn)行穩(wěn)定性的影響。研究結(jié)果表明：組合擋板能很好地降低液體對罐體結(jié)構(gòu)的沖擊強(qiáng)度。雖然剛性擋板能有效地抑制容器內(nèi)液體的晃動(dòng)，但是剛性擋板會(huì)增大罐車結(jié)構(gòu)的質(zhì)量，還會(huì)使罐內(nèi)的清理工作變得麻煩。隨著復(fù)合材料的不斷發(fā)展，人們發(fā)現(xiàn)彈性材料也可以制作成擋板。不僅可以起到抑制罐內(nèi)液體晃動(dòng)的作用，而且還可以減輕罐車結(jié)構(gòu)質(zhì)量。Duan 等[6]采用流固耦合方法分析了在加速過程中的液體晃動(dòng)現(xiàn)象、擋板的變形及受力情況，研究了在充液比為0.5 和0.8 時(shí)水平圓柱三維模型的液相分布。研究表明：充液比為0.8 的儲(chǔ)罐液體分布更穩(wěn)定。Zhang等[7]運(yùn)用SPH （Smoothed Particle Hydrodynamics） 法研究了二維矩形罐多彈性擋板及其組合擋板對液體晃動(dòng)的影響。研究表明：垂直或T 形擋板與水平擋板的組合在抑制晃動(dòng)方面要優(yōu)于單個(gè)擋板。Saghi 等[8]提出了一種水彈性模型，研究由具柔性壁面的梯形和矩形儲(chǔ)罐上的搖擺引起的晃動(dòng)載荷。結(jié)果表明：梯形儲(chǔ)罐比矩形儲(chǔ)罐有更好的抑制液體晃動(dòng)的性能。Yu 等[9]采用有限元分析和任意拉格朗日-歐拉（ALE）方法研究了擋板滲透率對液體晃動(dòng)的影響。研究發(fā)現(xiàn)一個(gè)最有效的抑制晃動(dòng)的臨界滲透系數(shù)，并且?guī)в锌锥磽醢宓膬?chǔ)罐內(nèi)會(huì)出現(xiàn)旋渦。Cho 等[10]基于速度勢的非線性有限元法研究在水平激勵(lì)下二維罐內(nèi)液體的大幅晃動(dòng)，研究表明：在罐體內(nèi)液體的晃動(dòng)量與防波板的設(shè)計(jì)參數(shù)有很強(qiáng)的相關(guān)性。Meng等[11]提出了一種半解析數(shù)學(xué)模型來研究矩形剛性容器與彈性擋板相互作用下液體的晃動(dòng)特性。研究表明薄擋板能提高液體晃動(dòng)的固有頻率。包文紅等[12]采用Fluent 軟件研究液罐車在不同加速度激勵(lì)下，液體晃動(dòng)對罐壁的沖擊力隨時(shí)間變化的規(guī)律。研究結(jié)果表明液罐車在減速階段擋板受到的沖擊力最大。徐曉美等[13]對液罐車在不同充液比和不同激勵(lì)時(shí)的穩(wěn)定性問題進(jìn)行了研究。研究結(jié)果顯示液罐車充液比為50%時(shí)，液體晃動(dòng)對車輛側(cè)傾穩(wěn)定性影響最大；而且隨著充液比的增大，罐體所受的力和力矩也隨之增大。Pozzetti 等[14]在VOF 算法的基礎(chǔ)上，提出了一種全新的三相流多尺度方法。該方法在考慮不同長度尺度的影響的同時(shí)還能大幅度減少計(jì)算負(fù)擔(dān)，研究結(jié)果顯示：三相流多尺度方法比流體體積法提供了更高的精度。

基于上述研究，本文在傳統(tǒng)剛性擋板的基礎(chǔ)上探索研究彈性擋板的阻尼特性，研究液罐車以橫向斜坡階躍制動(dòng)時(shí)罐內(nèi)液體晃動(dòng)對擋板產(chǎn)生的沖擊，探討液罐車在這種沖擊下彈性擋板的相對長度（H/R）和厚度對液體晃動(dòng)的影響。此項(xiàng)研究將對罐車內(nèi)液體晃動(dòng)的抑制研究提供一種新的思路。

1 計(jì)算模型

1.1 流體域的控制方程

本文使用開源軟件OpenFOAM（Open Field Operation And Manipulation）對罐體內(nèi)的液體晃動(dòng)現(xiàn)象進(jìn)行模擬。采用兩相流模型對液體晃動(dòng)過程中出現(xiàn)的氣相（氣體）和液相（液體）進(jìn)行模擬。由于涉及液體晃動(dòng)的液體流速很小，相應(yīng)的雷諾數(shù)很小，晃動(dòng)模型采用層流模型即可。涉及液體晃動(dòng)的控制方程為不可壓縮流體的連續(xù)性方程和Navier-Stokes動(dòng)量方程為：

式中：ρ為罐體內(nèi)液體密度；U為速度矢量；t為時(shí)間；p為壓力；τ為剪切應(yīng)力。

在兩相流模型中，涉及到捕捉氣液交界面的問題，本研究采用VOF （Volume of Fluid） 方法來捕捉自由液面。VOF 法引入F（x，y，t）函數(shù)，主要是跟蹤自由面的流體體積分?jǐn)?shù)，表示流體占所在網(wǎng)格體積的百分?jǐn)?shù)。當(dāng)F=0 時(shí)，網(wǎng)格內(nèi)全為氣體域；當(dāng)F=1 時(shí)，網(wǎng)格全處于液體域；當(dāng)0＜F＜1 時(shí)，表示網(wǎng)格一定含有自由液面。因此，函數(shù)F的求解，就可以跟蹤自由液面的位置。函數(shù)F隨時(shí)間變化的控制方程為：

式中：θ為部分單元體參數(shù)，其值在0 到1 之間；u、v分別為x、y方向上的流體速度分量。

由于圓柱壁面邊界處流體和固體沒有相對運(yùn)動(dòng)，所以采用無滑移邊界條件；在運(yùn)動(dòng)中通量修正保證彈性擋板面通量為0，所以采用運(yùn)動(dòng)的壁面邊界條件。

1.2 彈性擋板的控制方程

由結(jié)構(gòu)力學(xué)，彈性擋板的控制方程可以表達(dá)為：

式中：ρ為彈性擋板密度；d為位移；f為體積力矢量；σ為柯西應(yīng)力張量。

使用定常剪切模量參數(shù)輸入超彈模型Neo-Hookean，其彈性應(yīng)變能函數(shù)表達(dá)式為：

式中：μ為剪切模量；I1為第一偏應(yīng)變常量；J為彈性變形梯度的行列式；D1為材料的不可壓縮參數(shù)，詳細(xì)定義為：

式中：γ為泊松比。

1.3 數(shù)值模型建立

選取液罐橫向截面，建立一個(gè)半徑R=1.015 m、厚度為0.05 m的二維圓形剛性罐體（罐體的厚度為一個(gè)網(wǎng)格尺寸大小）。如圖1 所示，在罐體中心處放置高為H、寬為L的彈性擋板，采用了無量綱的方法設(shè)置了8 種彈性擋板的相對長度，其具體參數(shù)設(shè)置如表1 所示，彈性擋板的具體材料參數(shù)如表2 所示，力矩的參考點(diǎn)在罐體的底部點(diǎn)O′，罐體充液比為50%。二維圓形模型的構(gòu)建、網(wǎng)格的劃分和流體域數(shù)值模擬計(jì)算均在OpenFOAM 軟件中進(jìn)行，其中網(wǎng)格是在blockMeshDict 里構(gòu)建的，固體域數(shù)值模擬計(jì)算在Deal.II軟件中進(jìn)行的。

表1 計(jì)算參數(shù)設(shè)置

表2 彈性擋板的材料設(shè)置

圖1 罐體和擋板結(jié)構(gòu)圖

1.4 晃動(dòng)響應(yīng)分析

在充液比為50%時(shí)，如果液罐車發(fā)生變道制動(dòng)、轉(zhuǎn)向以及換道運(yùn)動(dòng)，其自由液面晃動(dòng)最為劇烈，更易發(fā)生側(cè)翻現(xiàn)象[15]。故本文選取罐體的充液比為50%。液罐車在轉(zhuǎn)向時(shí)的離心加速度用斜坡階躍的加速度進(jìn)行模擬，加速度隨時(shí)間變化的關(guān)系如式（7）所示。

式中：a(t)為加速度激勵(lì)；k為加速度斜坡的上升率；a是穩(wěn)態(tài)加速度的幅值，由于模擬的是液罐車的轉(zhuǎn)向運(yùn)動(dòng)，幅值取0.25g；β為圓弧過度參數(shù)，一般取0.2。加速度曲線如圖2所示。

圖2 圓滑斜坡階躍加速度激勵(lì)

液體載荷轉(zhuǎn)移量是通過動(dòng)態(tài)液體質(zhì)心（cg）與靜止?fàn)顟B(tài)下液體質(zhì)心之間的變化來評估的，其相應(yīng)的計(jì)算公式為：

式中：xcg和ycg分別為x軸和y軸的液體載荷轉(zhuǎn)移分量；x0和y0為靜止?fàn)顟B(tài)下液體質(zhì)心坐標(biāo)；V為液體總體積；Va是液體單元a的體積，xa和ya是單元a的質(zhì)心坐標(biāo)；ψ為液體域，如圖1所示。

對于任意單元b，用單元中心處壓力與單元邊界面面積的乘積來計(jì)算任意單元b對壁面的作用力，將該作用力在整個(gè)濕周（液體與壁面接觸的部分）進(jìn)行積分來計(jì)算作用在壁面的總晃動(dòng)力，其計(jì)算公式為：

式中：F為液體晃動(dòng)對邊界面的總作用力；p為任意單元b的壓力；s為單元與邊界接觸的面積。

根據(jù)力矩的定義，任意單元b相對于參考點(diǎn)O′的矢徑與單元b產(chǎn)生的作用力的叉乘即為該單元對參考點(diǎn)O′的力矩。對該力矩在整個(gè)濕周進(jìn)行積分即為晃動(dòng)液體作用在壁面的總晃動(dòng)力矩。相應(yīng)公式為：

式中：M為作用在罐體壁面的總力矩；l為單元到參考點(diǎn)O′的位置矢量。

2 網(wǎng)格無關(guān)性驗(yàn)證

一般來說網(wǎng)格劃分越細(xì)密，其仿真的結(jié)果越精確。但是實(shí)際上當(dāng)網(wǎng)格加密到一定的程度時(shí)，仿真的結(jié)果變化就很小，而且會(huì)導(dǎo)致仿真時(shí)間成本大幅增加。故為了綜合效率和精度一般認(rèn)為，在滿足庫朗數(shù)限制條件下，當(dāng)仿真數(shù)據(jù)隨網(wǎng)格尺寸的變小而只產(chǎn)生微小變化時(shí)，就達(dá)到網(wǎng)格無關(guān)性了。網(wǎng)格尺寸對液體晃動(dòng)仿真結(jié)果影響較大，故在仿真中選取了3種不同數(shù)量的網(wǎng)格。

針對充液比為50%的流體域，在OpenFOAM 中使用blockMeshDict 對模型進(jìn)行分塊劃分及O 型塊處理，每種網(wǎng)格的模擬是在ax= 0.25g的圓滑過渡的斜坡階躍加速度下進(jìn)行的，以模擬罐車轉(zhuǎn)向時(shí)的轉(zhuǎn)向加速度。設(shè)置了3 種不同數(shù)量的網(wǎng)格（31 528、44 038 和61 946），分別代表不同精細(xì)程度的網(wǎng)格。如圖3 中3 種網(wǎng)格之間的計(jì)算結(jié)果相差微小，為了較小的仿真時(shí)間選取網(wǎng)格數(shù)為31 528為計(jì)算網(wǎng)格。

圖3 網(wǎng)格無關(guān)性驗(yàn)證

3 二維橫向模型晃動(dòng)分析的結(jié)果和討論

3.1 彈性擋板形變對液體晃動(dòng)的影響

圖4 為彈性擋板中心處水平位移的最大形變量隨H/R的值和擋板厚度變化的趨勢圖。由圖可知，彈性擋板中心處的水平位移峰值隨著H/R的比值的增大而增大。從圖4 中還可以看出，隨著彈性擋板厚度的增加，其中心處對應(yīng)的水平位移峰值減小；并且在同一H/R值時(shí)，彈性擋板的厚度增加5 mm 時(shí)，其水平位移峰值下降約10 mm。原因是當(dāng)楊氏模量相同的時(shí)候，剛度隨彈性板厚度的增大而逐漸增大，彈性擋板抵抗變形的能力相應(yīng)增強(qiáng)。

圖4 彈性擋板水平位移峰值

為研究彈性擋板中心處的水平位移量和壓力變化與液體晃動(dòng)的關(guān)系。圖5 為10 mm 厚度的彈性擋板在板的相對長度（H/R）值分別為1.1、1.3、1.5、1.7 時(shí)，擋板中心處的水平位移量和壓力的時(shí)間歷程曲線。由圖可知，彈性擋板中心處位移隨時(shí)間歷程曲線與正弦曲線相似，這表明在液體沖擊作用下，彈性擋板在其平衡位置來回振動(dòng)，并且其振動(dòng)的幅值隨著H/R值的增大而增大。

從圖5 中還可以看出，彈性擋板中心處的壓力峰值主要集中在0.7 s 附近，即液體晃動(dòng)的第一個(gè)周期內(nèi)，這表明在斜坡階躍加速度作用下，液體晃動(dòng)產(chǎn)生的沖擊載荷在第一個(gè)周期內(nèi)最大，此后沖擊載荷逐漸減小。對比不同長度的彈性擋板下液體沖擊載荷的峰值可以發(fā)現(xiàn)彈性擋板的長度對其壓力峰值的影響較小。另外，觀察壓力的時(shí)間歷程曲線可以看出，“雙峰”現(xiàn)象較為明顯，表明液體晃動(dòng)運(yùn)動(dòng)出現(xiàn)了非線性現(xiàn)象。由于液體與彈性擋板的相互作用，晃動(dòng)液體在較短的時(shí)間內(nèi)多次與彈性擋板發(fā)生碰撞，彈性板的形變吸收了一部分液體晃動(dòng)動(dòng)能，使得壓強(qiáng)峰值分化，避免出現(xiàn)過大的液體沖擊力。

3.2 彈性擋板對液體載荷轉(zhuǎn)移量的影響

液體載荷轉(zhuǎn)移量直接影響作用在罐車上整體的力和力矩大小，是評價(jià)液罐車運(yùn)行穩(wěn)定性的重要指標(biāo)。仿真中罐體的激勵(lì)為ax= 0.25g，充液比為50%，設(shè)置了3 種厚度彈性擋板，即10 mm，15 mm，20 mm。彈性擋板的相對長度（H/R）分別設(shè)置為1.0、1.1、1.2、1.3、1.4、1.5、1.6、1.7 共8 個(gè)值。為研究彈性擋板對液體載荷轉(zhuǎn)移量的影響，表3 列出了罐體內(nèi)無擋板和設(shè)置不同彈性擋板的情況下，其液體載荷轉(zhuǎn)移量峰值及其下降率。

表3 不同彈性擋板液體載荷轉(zhuǎn)移量的峰值及與無擋板橫縱載荷轉(zhuǎn)移量峰值相對應(yīng)的下降率

無擋板的橫縱載荷轉(zhuǎn)移量峰值分別為x0=0.193 m、y0=0.046 m。根據(jù)表3顯示，相比較于無擋板的罐體，帶彈性擋板的罐體的橫向載荷轉(zhuǎn)移量峰值降低19.9%～25.9%，并且其相應(yīng)的縱向載荷轉(zhuǎn)移量峰值降低32.4%～45.2%。從表3還可以看出：H/R值從1.0增大到1.6時(shí)，隨著彈性擋板厚度的增大，其對應(yīng)縱向和橫向載荷轉(zhuǎn)移量峰值的下降率呈減小趨勢；當(dāng)H/R值增大到1.7時(shí)，情況正好相反，但載荷轉(zhuǎn)移量的變化幅度不大，在0.6%左右。該結(jié)果表明：在選擇彈性擋板作為防晃裝置時(shí)，減小擋板的厚度可以提高其對應(yīng)載荷轉(zhuǎn)移量峰值的下降率，即提高擋板抑制液體晃動(dòng)的效果。

3.3 力和力矩

圖6 為橫向加速度激勵(lì)下，作用在罐體上力的峰值隨H/R值的變化趨勢圖。由圖可知，和無擋板模型相比，彈性擋板可以很好地減少圓形罐體橫向晃動(dòng)力的峰值。原因是當(dāng)水波撞擊彈性擋板時(shí)，彈性擋板可以很好地把液體的動(dòng)能轉(zhuǎn)化為彈性勢能，當(dāng)彈性勢能積蓄到一定的時(shí)候，會(huì)與后續(xù)的水波發(fā)生對沖消耗。并且在圖5 的“雙峰”現(xiàn)象也顯示了彈性擋板在降低液體對罐體沖擊力上有很好的效果。

圖6 橫向晃動(dòng)力峰值

圖7 為罐體內(nèi)無擋板和放置彈性擋板時(shí)，液體晃動(dòng)產(chǎn)生的側(cè)傾力矩峰值隨板的相對長度的變化趨勢圖。對于放置彈性擋板的罐體，罐體內(nèi)液體晃動(dòng)產(chǎn)生的側(cè)傾力矩的峰值明顯的小于罐體內(nèi)無擋板時(shí)對應(yīng)的側(cè)傾力矩峰值，帶彈性擋板的罐車具有更好的運(yùn)行穩(wěn)定性。

圖7 側(cè)傾力矩峰值

從圖6～7 還可以看出，彈性擋板的厚度變化對橫向晃動(dòng)力和側(cè)傾力矩的峰值影響不大。故在后續(xù)使用彈性擋板作為阻尼裝置時(shí)，建議選取較小的厚度，這樣既可以在性能上滿足需求還可以實(shí)現(xiàn)擋板的輕量化。

3.4 罐體內(nèi)的速度矢量和壓強(qiáng)分析

圖8 是H/R值為1.0、擋板厚度為10 mm 時(shí)，不同時(shí)刻流體域的速度矢量分布和壓強(qiáng)分布圖。從圖中可以看出，在自由液面附近出現(xiàn)了破碎波的現(xiàn)象，由于液體與氣體之間存在粘性摩擦力，液體表面速度越快，與氣體之間的耗散就越大，其動(dòng)能的損失轉(zhuǎn)換為內(nèi)能，并且不增加壓力。即液體與氣體邊界交換的耗散會(huì)影響流體的勢能。

圖8 液體表面破碎的速度矢量分布和壓強(qiáng)分布（H/R=1，板厚10 mm）

圖9 為帶彈性擋板的圓形罐體內(nèi)的流體域的速度矢量分布和壓強(qiáng)分布圖，其中時(shí)間為1.36 s，板厚10 mm，板的相對長度分別為1.0、1.4和1.7。在圖5中1～2 s內(nèi)壓強(qiáng)的峰值隨著彈性擋板的相對長度的增加，壓強(qiáng)卻漸漸減小并趨近于零。此數(shù)據(jù)可以由圖9 中的現(xiàn)象解釋，從圖中可以發(fā)現(xiàn)：隨著彈性擋板相對長度的增加，其兩側(cè)液面的高度差逐漸減小，此時(shí)擋板中心處的左右壓強(qiáng)差也在減小，此外，在擋板右下角處出現(xiàn)了渦流，并且隨著彈性擋板相對長度的增加，渦流明顯增大，底部的壓強(qiáng)也隨之增加。如圖8～9 中液體產(chǎn)生的渦流會(huì)使內(nèi)摩擦力做功，造成進(jìn)一步的耗散。

圖9 不同彈性擋板的速度矢量分布和壓強(qiáng)分布（t=1.36 s，板厚10 mm）

4 結(jié)束語

本文在傳統(tǒng)剛性擋板的基礎(chǔ)上探索研究了彈性擋板對在斜坡階躍加速度激勵(lì)的作用下圓形罐體內(nèi)對液體晃動(dòng)響應(yīng)的影響，得到以下結(jié)論。

（1）整個(gè)物理模型中，從彈性擋板到罐體的耗散主要由彈性擋板的耗散、液體產(chǎn)生渦流時(shí)的耗散以及液體表面破碎波的阻尼3個(gè)區(qū)域組成。

（2）彈性擋板能有效地降低液體載荷轉(zhuǎn)移量的峰值，從而減小由液體晃動(dòng)作用在罐車上的側(cè)傾力矩，提高了車輛轉(zhuǎn)向時(shí)的穩(wěn)定性。

（3）較小的彈性擋板厚度對晃動(dòng)的抑制效果更佳。

（4）對于不同厚度的彈性擋板，其中心處的水平位移峰值都隨著其相對長度（H/R）的增大而呈現(xiàn)出增大趨勢；當(dāng)H/R的值固定時(shí)，隨著擋板厚度的增加，其中心處的水平位移峰值減小。