徐 凱,王瓊瑤,林國勉
(五邑大學智能制造學部,廣東 江門 529000)
液體晃動現(xiàn)象常出現(xiàn)在部分充液容器中,并且涉及到工程領域的各個方面,包括航空航天領域、陸運領域、海洋工程領域以及儲罐等。液體的劇烈晃動可能導致充液容器的結構破壞、液罐車側翻事故的頻發(fā)以及航天器發(fā)射的失敗等。因此,研究如何抑制充液容器內(nèi)液體的晃動對工程應用具有十分重要的意義。
自20 世紀60年代開始,Abramson 等[1]比較全面地開展了常重力環(huán)境下各種形狀的擋板的液體晃動實驗研究。Hasheminejad 和Mohammadi[2]使用線性勢理論的二維動力學方法研究了液體在無黏性和不可壓縮的條件下,水平圓柱擋板容器的3 種擋板安置方式對液體晃動的影響。結果表明:容器在高充液比情況時,底部安裝垂直擋板不是抑制部分填充容器內(nèi)液體晃動的最有效的手段。Yu和Chu[3]研究了部分充滿液罐車的液體晃動的非線性運動,研究結果表明:在充液比在90%時,瞬時側向力最小。Yan 和Rakheja[4]研究了部分充液的液罐車在有、無擋板的情況下其直線制動性能特征。研究結果顯示當罐體內(nèi)有擋板時,液罐車的制動性能明顯提高。Wang 等[5]基于VOF (Volume of Fluid) 模型研究了無擋板,傳統(tǒng)剛性擋板和組合擋板在橫、縱方向上對液罐車運行穩(wěn)定性的影響。研究結果表明:組合擋板能很好地降低液體對罐體結構的沖擊強度。雖然剛性擋板能有效地抑制容器內(nèi)液體的晃動,但是剛性擋板會增大罐車結構的質量,還會使罐內(nèi)的清理工作變得麻煩。隨著復合材料的不斷發(fā)展,人們發(fā)現(xiàn)彈性材料也可以制作成擋板。不僅可以起到抑制罐內(nèi)液體晃動的作用,而且還可以減輕罐車結構質量。Duan 等[6]采用流固耦合方法分析了在加速過程中的液體晃動現(xiàn)象、擋板的變形及受力情況,研究了在充液比為0.5 和0.8 時水平圓柱三維模型的液相分布。研究表明:充液比為0.8 的儲罐液體分布更穩(wěn)定。Zhang等[7]運用SPH (Smoothed Particle Hydrodynamics) 法研究了二維矩形罐多彈性擋板及其組合擋板對液體晃動的影響。研究表明:垂直或T 形擋板與水平擋板的組合在抑制晃動方面要優(yōu)于單個擋板。Saghi 等[8]提出了一種水彈性模型,研究由具柔性壁面的梯形和矩形儲罐上的搖擺引起的晃動載荷。結果表明:梯形儲罐比矩形儲罐有更好的抑制液體晃動的性能。Yu 等[9]采用有限元分析和任意拉格朗日-歐拉(ALE)方法研究了擋板滲透率對液體晃動的影響。研究發(fā)現(xiàn)一個最有效的抑制晃動的臨界滲透系數(shù),并且?guī)в锌锥磽醢宓膬迌?nèi)會出現(xiàn)旋渦。Cho 等[10]基于速度勢的非線性有限元法研究在水平激勵下二維罐內(nèi)液體的大幅晃動,研究表明:在罐體內(nèi)液體的晃動量與防波板的設計參數(shù)有很強的相關性。Meng等[11]提出了一種半解析數(shù)學模型來研究矩形剛性容器與彈性擋板相互作用下液體的晃動特性。研究表明薄擋板能提高液體晃動的固有頻率。包文紅等[12]采用Fluent 軟件研究液罐車在不同加速度激勵下,液體晃動對罐壁的沖擊力隨時間變化的規(guī)律。研究結果表明液罐車在減速階段擋板受到的沖擊力最大。徐曉美等[13]對液罐車在不同充液比和不同激勵時的穩(wěn)定性問題進行了研究。研究結果顯示液罐車充液比為50%時,液體晃動對車輛側傾穩(wěn)定性影響最大;而且隨著充液比的增大,罐體所受的力和力矩也隨之增大。Pozzetti 等[14]在VOF 算法的基礎上,提出了一種全新的三相流多尺度方法。該方法在考慮不同長度尺度的影響的同時還能大幅度減少計算負擔,研究結果顯示:三相流多尺度方法比流體體積法提供了更高的精度。
基于上述研究,本文在傳統(tǒng)剛性擋板的基礎上探索研究彈性擋板的阻尼特性,研究液罐車以橫向斜坡階躍制動時罐內(nèi)液體晃動對擋板產(chǎn)生的沖擊,探討液罐車在這種沖擊下彈性擋板的相對長度(H/R)和厚度對液體晃動的影響。此項研究將對罐車內(nèi)液體晃動的抑制研究提供一種新的思路。
本文使用開源軟件OpenFOAM(Open Field Operation And Manipulation)對罐體內(nèi)的液體晃動現(xiàn)象進行模擬。采用兩相流模型對液體晃動過程中出現(xiàn)的氣相(氣體)和液相(液體)進行模擬。由于涉及液體晃動的液體流速很小,相應的雷諾數(shù)很小,晃動模型采用層流模型即可。涉及液體晃動的控制方程為不可壓縮流體的連續(xù)性方程和Navier-Stokes動量方程為:
式中:ρ為罐體內(nèi)液體密度;U為速度矢量;t為時間;p為壓力;τ為剪切應力。
在兩相流模型中,涉及到捕捉氣液交界面的問題,本研究采用VOF (Volume of Fluid) 方法來捕捉自由液面。VOF 法引入F(x,y,t)函數(shù),主要是跟蹤自由面的流體體積分數(shù),表示流體占所在網(wǎng)格體積的百分數(shù)。當F=0 時,網(wǎng)格內(nèi)全為氣體域;當F=1 時,網(wǎng)格全處于液體域;當0<F<1 時,表示網(wǎng)格一定含有自由液面。因此,函數(shù)F的求解,就可以跟蹤自由液面的位置。函數(shù)F隨時間變化的控制方程為:
式中:θ為部分單元體參數(shù),其值在0 到1 之間;u、v分別為x、y方向上的流體速度分量。
由于圓柱壁面邊界處流體和固體沒有相對運動,所以采用無滑移邊界條件;在運動中通量修正保證彈性擋板面通量為0,所以采用運動的壁面邊界條件。
由結構力學,彈性擋板的控制方程可以表達為:
式中:ρ為彈性擋板密度;d為位移;f為體積力矢量;σ為柯西應力張量。
使用定常剪切模量參數(shù)輸入超彈模型Neo-Hookean,其彈性應變能函數(shù)表達式為:
式中:μ為剪切模量;I1為第一偏應變常量;J為彈性變形梯度的行列式;D1為材料的不可壓縮參數(shù),詳細定義為:
式中:γ為泊松比。
選取液罐橫向截面,建立一個半徑R=1.015 m、厚度為0.05 m的二維圓形剛性罐體(罐體的厚度為一個網(wǎng)格尺寸大?。?。如圖1 所示,在罐體中心處放置高為H、寬為L的彈性擋板,采用了無量綱的方法設置了8 種彈性擋板的相對長度,其具體參數(shù)設置如表1 所示,彈性擋板的具體材料參數(shù)如表2 所示,力矩的參考點在罐體的底部點O′,罐體充液比為50%。二維圓形模型的構建、網(wǎng)格的劃分和流體域數(shù)值模擬計算均在OpenFOAM 軟件中進行,其中網(wǎng)格是在blockMeshDict 里構建的,固體域數(shù)值模擬計算在Deal.II軟件中進行的。
表1 計算參數(shù)設置
表2 彈性擋板的材料設置
圖1 罐體和擋板結構圖
在充液比為50%時,如果液罐車發(fā)生變道制動、轉向以及換道運動,其自由液面晃動最為劇烈,更易發(fā)生側翻現(xiàn)象[15]。故本文選取罐體的充液比為50%。液罐車在轉向時的離心加速度用斜坡階躍的加速度進行模擬,加速度隨時間變化的關系如式(7)所示。
式中:a(t)為加速度激勵;k為加速度斜坡的上升率;a是穩(wěn)態(tài)加速度的幅值,由于模擬的是液罐車的轉向運動,幅值取0.25g;β為圓弧過度參數(shù),一般取0.2。加速度曲線如圖2所示。
圖2 圓滑斜坡階躍加速度激勵
液體載荷轉移量是通過動態(tài)液體質心(cg)與靜止狀態(tài)下液體質心之間的變化來評估的,其相應的計算公式為:
式中:xcg和ycg分別為x軸和y軸的液體載荷轉移分量;x0和y0為靜止狀態(tài)下液體質心坐標;V為液體總體積;Va是液體單元a的體積,xa和ya是單元a的質心坐標;ψ為液體域,如圖1所示。
對于任意單元b,用單元中心處壓力與單元邊界面面積的乘積來計算任意單元b對壁面的作用力,將該作用力在整個濕周(液體與壁面接觸的部分)進行積分來計算作用在壁面的總晃動力,其計算公式為:
式中:F為液體晃動對邊界面的總作用力;p為任意單元b的壓力;s為單元與邊界接觸的面積。
根據(jù)力矩的定義,任意單元b相對于參考點O′的矢徑與單元b產(chǎn)生的作用力的叉乘即為該單元對參考點O′的力矩。對該力矩在整個濕周進行積分即為晃動液體作用在壁面的總晃動力矩。相應公式為:
式中:M為作用在罐體壁面的總力矩;l為單元到參考點O′的位置矢量。
一般來說網(wǎng)格劃分越細密,其仿真的結果越精確。但是實際上當網(wǎng)格加密到一定的程度時,仿真的結果變化就很小,而且會導致仿真時間成本大幅增加。故為了綜合效率和精度一般認為,在滿足庫朗數(shù)限制條件下,當仿真數(shù)據(jù)隨網(wǎng)格尺寸的變小而只產(chǎn)生微小變化時,就達到網(wǎng)格無關性了。網(wǎng)格尺寸對液體晃動仿真結果影響較大,故在仿真中選取了3種不同數(shù)量的網(wǎng)格。
針對充液比為50%的流體域,在OpenFOAM 中使用blockMeshDict 對模型進行分塊劃分及O 型塊處理,每種網(wǎng)格的模擬是在ax= 0.25g的圓滑過渡的斜坡階躍加速度下進行的,以模擬罐車轉向時的轉向加速度。設置了3 種不同數(shù)量的網(wǎng)格(31 528、44 038 和61 946),分別代表不同精細程度的網(wǎng)格。如圖3 中3 種網(wǎng)格之間的計算結果相差微小,為了較小的仿真時間選取網(wǎng)格數(shù)為31 528為計算網(wǎng)格。
圖3 網(wǎng)格無關性驗證
圖4 為彈性擋板中心處水平位移的最大形變量隨H/R的值和擋板厚度變化的趨勢圖。由圖可知,彈性擋板中心處的水平位移峰值隨著H/R的比值的增大而增大。從圖4 中還可以看出,隨著彈性擋板厚度的增加,其中心處對應的水平位移峰值減??;并且在同一H/R值時,彈性擋板的厚度增加5 mm 時,其水平位移峰值下降約10 mm。原因是當楊氏模量相同的時候,剛度隨彈性板厚度的增大而逐漸增大,彈性擋板抵抗變形的能力相應增強。
圖4 彈性擋板水平位移峰值
為研究彈性擋板中心處的水平位移量和壓力變化與液體晃動的關系。圖5 為10 mm 厚度的彈性擋板在板的相對長度(H/R)值分別為1.1、1.3、1.5、1.7 時,擋板中心處的水平位移量和壓力的時間歷程曲線。由圖可知,彈性擋板中心處位移隨時間歷程曲線與正弦曲線相似,這表明在液體沖擊作用下,彈性擋板在其平衡位置來回振動,并且其振動的幅值隨著H/R值的增大而增大。
從圖5 中還可以看出,彈性擋板中心處的壓力峰值主要集中在0.7 s 附近,即液體晃動的第一個周期內(nèi),這表明在斜坡階躍加速度作用下,液體晃動產(chǎn)生的沖擊載荷在第一個周期內(nèi)最大,此后沖擊載荷逐漸減小。對比不同長度的彈性擋板下液體沖擊載荷的峰值可以發(fā)現(xiàn)彈性擋板的長度對其壓力峰值的影響較小。另外,觀察壓力的時間歷程曲線可以看出,“雙峰”現(xiàn)象較為明顯,表明液體晃動運動出現(xiàn)了非線性現(xiàn)象。由于液體與彈性擋板的相互作用,晃動液體在較短的時間內(nèi)多次與彈性擋板發(fā)生碰撞,彈性板的形變吸收了一部分液體晃動動能,使得壓強峰值分化,避免出現(xiàn)過大的液體沖擊力。
液體載荷轉移量直接影響作用在罐車上整體的力和力矩大小,是評價液罐車運行穩(wěn)定性的重要指標。仿真中罐體的激勵為ax= 0.25g,充液比為50%,設置了3 種厚度彈性擋板,即10 mm,15 mm,20 mm。彈性擋板的相對長度(H/R)分別設置為1.0、1.1、1.2、1.3、1.4、1.5、1.6、1.7 共8 個值。為研究彈性擋板對液體載荷轉移量的影響,表3 列出了罐體內(nèi)無擋板和設置不同彈性擋板的情況下,其液體載荷轉移量峰值及其下降率。
表3 不同彈性擋板液體載荷轉移量的峰值及與無擋板橫縱載荷轉移量峰值相對應的下降率
無擋板的橫縱載荷轉移量峰值分別為x0=0.193 m、y0=0.046 m。根據(jù)表3顯示,相比較于無擋板的罐體,帶彈性擋板的罐體的橫向載荷轉移量峰值降低19.9%~25.9%,并且其相應的縱向載荷轉移量峰值降低32.4%~45.2%。從表3還可以看出:H/R值從1.0增大到1.6時,隨著彈性擋板厚度的增大,其對應縱向和橫向載荷轉移量峰值的下降率呈減小趨勢;當H/R值增大到1.7時,情況正好相反,但載荷轉移量的變化幅度不大,在0.6%左右。該結果表明:在選擇彈性擋板作為防晃裝置時,減小擋板的厚度可以提高其對應載荷轉移量峰值的下降率,即提高擋板抑制液體晃動的效果。
圖6 為橫向加速度激勵下,作用在罐體上力的峰值隨H/R值的變化趨勢圖。由圖可知,和無擋板模型相比,彈性擋板可以很好地減少圓形罐體橫向晃動力的峰值。原因是當水波撞擊彈性擋板時,彈性擋板可以很好地把液體的動能轉化為彈性勢能,當彈性勢能積蓄到一定的時候,會與后續(xù)的水波發(fā)生對沖消耗。并且在圖5 的“雙峰”現(xiàn)象也顯示了彈性擋板在降低液體對罐體沖擊力上有很好的效果。
圖6 橫向晃動力峰值
圖7 為罐體內(nèi)無擋板和放置彈性擋板時,液體晃動產(chǎn)生的側傾力矩峰值隨板的相對長度的變化趨勢圖。對于放置彈性擋板的罐體,罐體內(nèi)液體晃動產(chǎn)生的側傾力矩的峰值明顯的小于罐體內(nèi)無擋板時對應的側傾力矩峰值,帶彈性擋板的罐車具有更好的運行穩(wěn)定性。
圖7 側傾力矩峰值
從圖6~7 還可以看出,彈性擋板的厚度變化對橫向晃動力和側傾力矩的峰值影響不大。故在后續(xù)使用彈性擋板作為阻尼裝置時,建議選取較小的厚度,這樣既可以在性能上滿足需求還可以實現(xiàn)擋板的輕量化。
圖8 是H/R值為1.0、擋板厚度為10 mm 時,不同時刻流體域的速度矢量分布和壓強分布圖。從圖中可以看出,在自由液面附近出現(xiàn)了破碎波的現(xiàn)象,由于液體與氣體之間存在粘性摩擦力,液體表面速度越快,與氣體之間的耗散就越大,其動能的損失轉換為內(nèi)能,并且不增加壓力。即液體與氣體邊界交換的耗散會影響流體的勢能。
圖8 液體表面破碎的速度矢量分布和壓強分布(H/R=1,板厚10 mm)
圖9 為帶彈性擋板的圓形罐體內(nèi)的流體域的速度矢量分布和壓強分布圖,其中時間為1.36 s,板厚10 mm,板的相對長度分別為1.0、1.4和1.7。在圖5中1~2 s內(nèi)壓強的峰值隨著彈性擋板的相對長度的增加,壓強卻漸漸減小并趨近于零。此數(shù)據(jù)可以由圖9 中的現(xiàn)象解釋,從圖中可以發(fā)現(xiàn):隨著彈性擋板相對長度的增加,其兩側液面的高度差逐漸減小,此時擋板中心處的左右壓強差也在減小,此外,在擋板右下角處出現(xiàn)了渦流,并且隨著彈性擋板相對長度的增加,渦流明顯增大,底部的壓強也隨之增加。如圖8~9 中液體產(chǎn)生的渦流會使內(nèi)摩擦力做功,造成進一步的耗散。
圖9 不同彈性擋板的速度矢量分布和壓強分布(t=1.36 s,板厚10 mm)
本文在傳統(tǒng)剛性擋板的基礎上探索研究了彈性擋板對在斜坡階躍加速度激勵的作用下圓形罐體內(nèi)對液體晃動響應的影響,得到以下結論。
(1)整個物理模型中,從彈性擋板到罐體的耗散主要由彈性擋板的耗散、液體產(chǎn)生渦流時的耗散以及液體表面破碎波的阻尼3個區(qū)域組成。
(2)彈性擋板能有效地降低液體載荷轉移量的峰值,從而減小由液體晃動作用在罐車上的側傾力矩,提高了車輛轉向時的穩(wěn)定性。
(3)較小的彈性擋板厚度對晃動的抑制效果更佳。
(4)對于不同厚度的彈性擋板,其中心處的水平位移峰值都隨著其相對長度(H/R)的增大而呈現(xiàn)出增大趨勢;當H/R的值固定時,隨著擋板厚度的增加,其中心處的水平位移峰值減小。