王瓊瑤， Subhash Rakheja， 上官文斌

(華南理工大學 機械與汽車工程學院， 廣州 510640)

貯箱內(nèi)液體晃動問題是工程中常見的問題，液罐車，水壩，航天器的燃料貯箱等都存在液體晃動的問題[1-2]。液罐車為液體貨物道路運輸?shù)淖畲筝d體，主要承擔石油、液化天然氣、危險化學品的運輸。由于液體貨物密度的不同以及車軸最大載荷的限制等原因，液罐車通常處于部分充液狀態(tài)。罐車轉(zhuǎn)向時，若行進速度較大，罐體內(nèi)的液體將大幅度晃動并與罐車相互作用，從而降低了液罐車的側(cè)傾穩(wěn)定性[3]。事故統(tǒng)計表明，相比較于固體載貨車輛，液罐車更容易發(fā)生側(cè)翻事故，且造成嚴重的人員傷亡和財產(chǎn)損失[4]。因此研究液罐車罐體內(nèi)液體晃動的動態(tài)特性以及如何限制液體晃動具有重要的現(xiàn)實意義。

限制罐體內(nèi)液體晃動的主要方法一般是在罐體內(nèi)添加隔板[5-7]或者防波板[8-10](帶孔的隔板)。Rakheja等采用了準靜態(tài)模型法研究了隔板的數(shù)目對液罐車側(cè)傾極限的影響。研究結(jié)果表明隔板能夠極大地提高液罐車的側(cè)傾穩(wěn)定性并減小罐體的部分充液狀態(tài)。Popov等首次基于全非線性流體動力學方程研究了帶隔板的罐體內(nèi)液體復雜的晃動現(xiàn)象。結(jié)果表明，與防波板相比，隔板能使罐車轉(zhuǎn)向時的側(cè)傾力矩減小30%。Wang等對隔板的研究表明，在罐車直線制動過程中，等間距的隔板能使液體貨物的縱向載荷轉(zhuǎn)移量最小。劉小民等采用雙向流固耦合的方法，對運動罐車在制動過程中罐體內(nèi)液體的晃動現(xiàn)象進行了數(shù)值分析，研究了不同材料的防波板在載液罐體內(nèi)的受力情況。研究表明，罐體端面受力隨著充液比的增大而增加，當充液比為0.8時，罐體的總體受力最大。Kandasamy等研究了傳統(tǒng)防波板、部分式防波板以及傾斜式防波板對液體晃動的影響。結(jié)果表明，傳統(tǒng)防波板在罐車制動過程中能有效抑制液體的縱向晃動，但在罐車轉(zhuǎn)向時對液體的橫向晃動影響不大。傾斜式防波板既能抑制液體的縱向晃動，也能抑制液體的橫向晃動。盡管隔板或者防波板能夠有效限制液體晃動，但它們極大地增大了罐車的結(jié)構(gòu)重量從而減小了運輸效率。此外，隔板或者防波板的引入使得罐體內(nèi)的清潔工作非常困難。

為了限制液體的晃動，一種新的抑制液體晃動的方法是在液體的自由面上覆蓋一層彈性膜，彈性膜的邊緣固定在罐體的內(nèi)壁。相比較于隔板或者防波板，彈性膜結(jié)構(gòu)簡單，重量輕，能有效提高液罐車的運輸效率。Bauer[11]研究了豎直圓柱形罐體內(nèi)液體晃動的共振頻率，其中液體表面覆蓋了一層彈性膜，且假設液體為不可壓縮無黏性的流體。研究結(jié)果表明，流體-彈性膜的耦合頻率大于不耦合時液體的晃動頻率。隨后，Bauer等[12]進一步研究了液體的黏性對流體-彈性膜的耦合頻率的影響，結(jié)果表明，液體的黏性會進一步增大耦合頻率。Parasil等[13-14]采用了任意拉格朗日歐拉方法研究了二維/三維矩形槽內(nèi)液體晃動的動力學行為，其中液體的自由面覆蓋了一層彈性膜。研究結(jié)果表明彈性膜本身的非線性特性對液體的晃動的影響非常大。

以上關于彈性膜對液體晃動的影響的研究僅限于豎直放置的罐體，且研究中僅考慮了液體與彈性膜的耦合作用，并沒有考慮空氣介質(zhì)。本文建立了流體-彈性膜耦合動力學模型(模型中考慮了空氣介質(zhì))，通過實驗驗證了模型的有效性。研究了二維水平圓柱形罐體在橫向加速度激勵下液體晃動的動力學現(xiàn)象。通過對比不帶彈性膜的罐體內(nèi)液體晃動的動力學響應，探究了彈性膜對液體晃動響應的影響。

1 流體-彈性膜耦合動力學模型

為了限制液體的晃動，本文提出了一種新型的抑制液體晃動的裝置：彈性膜，即在自由液面上覆蓋一層彈性膜，彈性膜的四周固定在罐體的內(nèi)壁。圖1為帶彈性膜的液灌模型圖。其中，罐體的橫截面為圓形，直徑為2.03 m。

圖1 液體二維橫向晃動的流體-彈性膜耦合模型

流體域采用的是二相流模型，即可壓縮的氣相和輕微可壓縮的液相，氣相和液相由彈性膜分開。在圓心位置建立笛卡爾坐標系，其中Y軸沿罐體的徑向方向，Z軸豎直向上。流體-彈性膜耦合問題是通過將計算流體動力學與流固耦合仿真平臺相結(jié)合來進行模擬的。在此流固耦合模型中，考慮了兩個流固耦合面，即彈性膜的上表面與氣相的交界面(流固耦合面1)以及彈性膜的下表面與液相的交界面(流固耦合面2)；與液相接觸的罐體部分采用壁面邊界條件；與氣相接觸的罐體部分則采用自由邊界條件，允許空氣自由出入罐體，如圖1所示。彈性膜的初始狀態(tài)假設為無預應力，且其固定于罐體內(nèi)壁的位置高度為當充液比(液體體積與罐體總體積的比值)為48%時自由液面的位置高度。因此，當充液體積比大于48%時，彈性膜將承受一定的預應力。

考慮到液體的黏性力與慣性力相比非常小，因此在建模中忽略了液體的黏性力。另外，由于液體晃動的過程中彈性膜較大的阻尼作用，導致液體晃動的幅值比較小，因此采用層流模型來模擬液體的晃動。彈性膜結(jié)構(gòu)在ANSYS中被劃分成結(jié)構(gòu)網(wǎng)格。流體域的網(wǎng)格劃分則通過CFX中的網(wǎng)格自動劃分技術完成。由于在彈性膜鄰近的網(wǎng)格變形量大，物理量的變化梯度大，因此對彈性膜鄰近的網(wǎng)格進行局部網(wǎng)格細化處理。由于是二維分析，垂直于紙面的方向的長度恰好為一個網(wǎng)格單元的長度。

通過對流固耦合模型的動力學響應進行分析，可以確定液體晃動產(chǎn)生的作用在罐體上的力的瞬態(tài)和穩(wěn)態(tài)值、液體載荷的轉(zhuǎn)移量以及側(cè)傾力矩的大小。液體載荷轉(zhuǎn)移量是用液體的質(zhì)心位置相對于液體處于靜態(tài)時質(zhì)心位置的變化量來表示的，如圖2所示。對流體域進行網(wǎng)格劃分(離散化)，則流體域液相的質(zhì)心坐標變化可通過流體域液相的體積積分求得

(1)

對于單元e，其對罐體壁面產(chǎn)生的作用力可由單元中心的壓強Pe與單元的面積矢量Ae求出。因而，液體晃動產(chǎn)生的合力為

(2)

式中：Fi為晃動合力沿i(i=y,z)軸的分量，Aei為單元e的面積矢量沿i軸的投影面積。φ為液相與罐體壁面以及彈性膜的接觸區(qū)域。

圖2 液體的載荷轉(zhuǎn)移量，力以及力矩的計算原理圖

側(cè)傾力矩的計算參考點為罐體幾何中心點在罐體底部的投影點O′,通過單個單元的力矩在區(qū)域φ上的積分求得

(3)

式中：Fe為單元e產(chǎn)生的晃動力矢量，le為單元相對于O′的位置矢量，M為側(cè)傾力矩。

2 模型有效性檢驗

圖3為帶彈性膜的罐體實驗模型，其中罐體長1.99 m，半徑0.268 m。罐體通過支撐架與支撐板固定在振動臺上。實驗用的液體為水，充液比為50%。

(a) 帶彈性氣袋的罐體模型實驗圖

(b) 帶彈性氣袋的罐體

在實驗中，由于氣密性的問題，故采用軟質(zhì)PVC材料制成的彈性氣袋來替代彈性膜，為了更好地模擬氣袋的密閉性，將與空氣介質(zhì)接觸的罐體定義為壁面邊界，而并非圖1所示的開放邊界。彈性氣袋的橫截面結(jié)構(gòu)如圖3(b)所示，整個彈性氣袋的形狀類似于一個半圓柱體(以滿足實驗中充液體積比為50%的設計要求)，長度和半徑均與罐體相同，以便于彈性氣袋固定在罐體的上圓柱面上。彈性氣袋通過尼龍搭扣被初步固定在罐體內(nèi)壁上，然后通過工業(yè)強力膠水對氣袋與罐體接觸的邊緣進行黏貼密封處理，其與液面接觸的部分因此可以近似看成是一層彈性膜。

為了驗證流體-彈性膜耦合動力學模型的有效性，首先需獲得實驗中彈性膜的材料參數(shù)。為此，對實驗用的制成彈性氣袋的材料樣本在MTS實驗臺上進行了多次拉伸實驗，如圖4(a)所示，MTS實驗臺的拉伸速率恒定為0.1 mm/s。將測試結(jié)果取平均值即得到圖4(b)的材料樣本的應力-應變曲線。

圖4 彈性膜材料拉伸實驗示意圖

從圖4中材料的應力-應變曲線關系可以看出，實驗用的材料可以看成是超彈性材料。在流固耦合模型中，采用Yeoh超彈性材料模型對彈性膜進行模擬。Yeoh模型主要用來模擬諸如橡膠等輕微可壓縮、非線性彈性材料的變形。模型的彈性屬性主要用應變能密度函數(shù)來表示，應變能密度函數(shù)一般用關于應變不變量I1，I2，I3的多項式表示。當模型用來模擬完全不可壓縮的橡膠材料時，相應的應變能密度函數(shù)僅是關于應變不變量I1的多項式。此時的應變能密度函數(shù)可以寫成

(4)

其中，應變不變量I1可以寫成

(5)

式中：Ci為材料常數(shù)，λ1，λ2，λ3是主伸長比。利用ANSYS WORKBENCH中的超彈性曲線擬合工具，將實驗得到的材料樣本的應力-應變曲線擬合成Yeoh超彈性材料模型。相應的材料常數(shù)便可確定：C1=2.310 8 Pa，C2=0.549 6 Pa，C3=0.120 8 Pa。

在實驗中，通過對液體橫向自由晃動響應進行分析可知，液體橫向晃動的固有頻率為1.3 Hz。為了驗證仿真模型的有效性，需對仿真模型在固有頻率處進行驗證。圖5為實驗結(jié)果與仿真結(jié)果的對比，其中橫向加速度激勵為ay=0.5sin(2πft)，激勵頻率為1.3 Hz，激勵的幅值分別為0.5 m/s2和1.0 m/s2。由于仿真模型是二維模型，因此實驗結(jié)果和仿真結(jié)果均是相對于單位罐體長度后的結(jié)果。由圖可知，仿真結(jié)果和實驗結(jié)果基本吻合，但在峰值處還是有一些偏差，主要有以下幾個原因：① 由于空間狹小，操作不便，彈性氣袋邊緣與罐體壁面處的黏貼情況在一些地方可能并不牢固；② 實驗中，液體晃動的過程中，彈性氣袋的氣壓會發(fā)生變化，而在仿真模型中，氣壓則不會變化。

3 流體-彈性膜耦合動力學模型的仿真分析

以實際液罐車的罐體尺寸為依據(jù)，建立了1∶1的簡化幾何模型。罐體的橫截面為圓形，直徑為2.03 m，由于是二維分析，垂直于紙面的方向的長度恰好為一個網(wǎng)格單元的長度：0.02 m。在液體的自由面上覆蓋一層彈性膜，彈性膜的兩端固定在罐體的內(nèi)壁上，彈性膜的厚度δ為0.005 m?？紤]了三種不同的充液體積比進行研究，即50%，65%和80%。液體載荷在不同的充液體積比時的密度并不相同，以保持液體載荷的總重量不變。不同充液體積比下，液體載荷的密度以及代表性的液體，如表1所示。

(a) 橫向加速度激勵：ay=0.5sin(2πft)，f=1.3 Hz

(b) 橫向加速度激勵：ay=1.0sin(2πft)，f=1.3 Hz

由于彈性膜的固定位置為充液高度為當充液比為48%時自由液面的位置高度，因此罐體內(nèi)的初始充液體積比為48%。這意味著在這三種不同充液體積比的情況下，彈性膜的變形大小將會不同，其預應力的大小也將不同，圖6為彈性膜在這三種充液體積比時的變形圖。

為了達到目標充液體積比(50%、65%和80%)，則需要往罐體內(nèi)額外注入液體。因此本文中一個完整的仿真模擬一般包含以下兩個階段：

(1) 充液階段：由于罐體內(nèi)的初始充液體積比為48%，為了達到仿真模擬所需的50%，65%和80%的充液體積比，需要往罐體內(nèi)注入額外的液體質(zhì)量Δm，Δm=ρΔV，ρ為不同充液體積比時液體的質(zhì)量密度；ΔV為仿真模擬所需的液體體積與罐體內(nèi)初始液體體積的差值。為了往罐體內(nèi)注入額外的液體質(zhì)量，在罐體接近底部的位置設置一個質(zhì)量點源，通過質(zhì)量點源流進罐體內(nèi)的質(zhì)量流量如圖7(a)所示，橫坐標表示時間(s)，縱坐標表示質(zhì)量流量(kg/s)，曲線與橫坐標圍成的面積即為往罐體額外注入的液體質(zhì)量Δm。

(2) 加速度激勵階段：當充液階段完成后，開始對罐體的流體域施加圓滑過渡的斜坡階躍加速度激勵，如圖7(b)所示，其數(shù)學表達式為

G(t)=

(6)

式中：G(t)表示加速度激勵，k表示上升速率，β表示圓滑過渡參數(shù)，一般取值0.2。

圖7 一次完整仿真的兩個階段

由于仿真中加速度激勵的大小為0.3g，液體晃動時彈性膜的應變不大，由圖4(b)可知。且本文提出的模型為二維模型，因此，在應變較小的情況下，彈性膜的楊氏模量可近似為恒定值，以提高計算效率。因此在本文中，彈性膜假設為由具有各向同性的橡膠材料制成，其物理參數(shù)如表2所示。

由于在仿真計算中，彈性膜的楊氏模量為一恒定值1×107Pa，因此在液體晃動的過程中，彈性膜的應力σ可以通過其應變ε來求得

σ=E·ε

(7)

在二維分析中，彈性膜的應變即為彈性膜的伸長率—彈性膜的伸長長度ΔL與其初始長度L0的比值

(8)

(i=1,2,…,N,N→∞)

(9)

彈性膜的總長可以由以下表達式計算

(10)

圖8 彈性膜的變形示意圖

3.1 晃動力和力矩

圖9分別為罐體內(nèi)在有無彈性膜的情況下，橫向晃動力的時間歷程響應的對比和側(cè)傾力矩的時間歷程響應的對比。其中，充液體積比為65%，橫向加速度激勵圓滑過渡的為斜坡-階躍加速度激勵，大小為0.3g，如圖7(b)所示。對比結(jié)果表明，罐體內(nèi)在有無彈性膜的情況下，橫向晃動力的峰值和穩(wěn)態(tài)值大小相當(圖9(a))；但彈性膜使作用在罐體上的側(cè)傾力矩的峰值和穩(wěn)態(tài)均值均大幅減小(圖9(b))，主要是因為彈性膜限制了液體自由面的運動，相應地便減小了液體的載荷轉(zhuǎn)移量；此外，通過對液體晃動的穩(wěn)態(tài)響應進行頻譜分析可知：在罐體內(nèi)有無彈性膜的情況下，液體橫向晃動的固有頻率分別為1.35 Hz和0.62 Hz，說明彈性膜明顯增大了液體橫向晃動的固有頻率。

(a) 橫向晃動力

(b) 側(cè)傾力矩

圖10和11分別為罐體內(nèi)在有無彈性膜的情況下，橫向力峰值和穩(wěn)態(tài)均值的對比以及側(cè)傾力矩峰值和穩(wěn)態(tài)均值的對比。其中，加速度激勵的幅值為0.3g，充液體積比分別為50%，65%和80%。由圖可知，彈性膜使作用在罐體上的橫向晃動力的峰值和穩(wěn)態(tài)均值均增大了，主要是因為彈性膜限制自由面向上運動，導致聚集的液體對罐體壁面的壓強增大。由圖10(a)可知，當充液體積比為50%時，彈性膜的初始應力較小，使得自由液面的運動幅值相對較大，自由液面較大幅值的運動反過來又會使彈性膜產(chǎn)生較大額外應力。當充液體積比為65%或80%時，由于彈性膜的初始應力較大，使得自由液面的運動幅值相對較小，自由液面的運動使彈性膜產(chǎn)生的額外應力較小。因此，在有彈性膜的罐體內(nèi)，充液體積比為50%時，作用在罐體上的橫向晃動力的峰值最大。此外，通過對比圖10(a)和(b)可知，當充液體積比為65%或80%時，由于彈性膜的初始應力較大，彈性膜產(chǎn)生的額外應力較小，使得橫向晃動力的峰值和穩(wěn)態(tài)均值之間的差異較小。

(a) 峰值

(b) 穩(wěn)態(tài)均值

(a) 峰值

(b) 穩(wěn)態(tài)均值

由圖11可知，彈性膜使得作用在罐體上的側(cè)傾力矩的峰值和穩(wěn)態(tài)均值均明顯減小。主要原因是側(cè)傾力矩受橫向和垂向載荷轉(zhuǎn)移量的影響很大，而彈性膜極大地限制了液體晃動時的載荷轉(zhuǎn)移量。以充液體積比為50%的情況為例，圖12為罐體內(nèi)有、無彈性膜的情況下，當充液體積比為50%時，橫向和縱向載荷轉(zhuǎn)移量的時間歷程曲線。由圖可知，在彈性膜的阻尼作用下，液體晃動時的橫向和縱向載荷轉(zhuǎn)移量均大幅減小。此外，從時間歷程響應曲線可知，彈性膜增大了液體晃動響應的振蕩頻率，縮短了液體晃動進入穩(wěn)態(tài)時所需的時間。

(a) 橫向載荷轉(zhuǎn)移量

(b) 垂向載荷轉(zhuǎn)移量

3.2 晃動頻率

對液體自由晃動響應進行頻率分辨率為0.01 Hz的快速傅里葉變換即可得到液體晃動的固有頻率。當充液體積比的范圍在50%～80%時，在帶彈性模的罐體內(nèi)，液體橫向晃動的固有頻率是不帶彈性膜罐體內(nèi)液體橫向晃動固有頻率的2倍至2.5倍，如圖13所示。當充液體積比從50%增大到80%時，不帶彈性膜罐體內(nèi)液體橫向晃動的固有頻率增加了15.5%，而帶彈性膜的罐體內(nèi)液體橫向晃動的固有頻率增加了30%，主要是高充液體積比下，彈性膜的應力值迅速增大引起的。由圖可知，在一般充液體積比下，帶彈性膜罐體內(nèi)液體橫向晃動的固有頻率在1.16～1.51 Hz范圍內(nèi)，而液罐車的緊急轉(zhuǎn)向頻率一般在0.6 Hz作用。因此，彈性膜使液體橫向晃動的固有頻率遠離了液罐車的緊急轉(zhuǎn)向頻率，從而避免了液罐車緊急轉(zhuǎn)向時共振現(xiàn)象的發(fā)生。

3.3 彈性膜的變形分析

圖14分別為彈性膜在液體晃動的初始狀態(tài)、穩(wěn)態(tài)以及彈性膜達到最大變形量時的輪廓圖。當充液體積比為50%時，由于彈性膜的應力較小，彈性膜的運動主要受液體運動的影響且彈性膜對液體晃動的影響相對

圖13 液體橫向晃動的固有頻率

較小。因此，在耦合模型中，除了晃動幅值相對較小外，液體的運動方式與晃動模型中液體的運動方式一致。值得注意的是，當晃動達到穩(wěn)態(tài)時，彈性膜關于其初始位置對稱。當充液體積比達到65%或者80%時，由于彈性膜的變形量急劇增大，其初始應力急劇增大，使得彈性膜在此變形基礎上的進一步變形量非常小，如圖14(b)和(c)所示。

(a) 充液體積比：50%

(b) 充液體積比：65%

(c) 充液體積比：80%

圖15為彈性膜各個狀態(tài)下的應力值隨充液體積比的變化趨勢。由圖可知，隨著充液體積比的增大，彈性膜的在這三種狀態(tài)下的應力值都急劇增大，以晃動初始狀態(tài)時彈性膜的應力值為例，當充液體積比為50%時，彈性膜的初始應力值為0.03 MPa，而當充液體積比為80%時，這個應力值增大到了2.95 MPa，增大了97倍。彈性膜在液體晃動的穩(wěn)態(tài)時、彈性膜變形量最大時的應力值隨著充液體積比的變化趨勢也與之類似。但是，隨著充液體積比的增大，彈性膜的應力值在這三種狀態(tài)之間的變化卻越來越?。寒敵湟后w積比為50%時，彈性膜的應力值從液體晃動初始狀態(tài)時的0.03 MPa到彈性膜變形量最大時的0.68 MPa，應力值的增長率達到了2 167%；當充液體積為65%和80%時，相應的應力值的增長率分別為19.2%和2.0%。

圖15 彈性膜在不通充液體積比、不同狀態(tài)下的應力值

4 結(jié) 論

(1) 建立了帶彈性膜罐體內(nèi)液體晃動的流體-彈性膜耦合動力學模型。并通過實驗驗證了模型的有效性。

(2) 研究了彈性膜及其預應力的大小對液體晃動響應的影響。研究結(jié)果表明，彈性膜能限制液體晃動的幅值，從而減小液體晃動對罐體壁面產(chǎn)生的側(cè)傾力矩。

(3) 彈性膜增大了液體橫向晃動的固有頻率。當充液體積比從50%增大到80%時，不帶彈性膜罐體內(nèi)液體橫向晃動的固有頻率增加了15.5%，而帶彈性膜的罐體內(nèi)液體橫向晃動的固有頻率增加了30%。此外，彈性膜使液體橫向晃動的固有頻率遠離了液罐車的緊急轉(zhuǎn)向頻率，從而避免了液罐車緊急轉(zhuǎn)向時共振現(xiàn)象的發(fā)生。