基于VOF法的球罐內(nèi)液體晃蕩數(shù)值模擬及載荷計算

劉岑凡 楊智榮 于哲敏 王國軍

(1.中國特種設(shè)備檢測研究院 北京 100029)(2.大連理工大學(xué) 大連 116024)

特種設(shè)備中許多壓力容器是液體的儲存裝置，如各種儲罐、球罐和槽車等。當發(fā)生地震或其他載荷激勵時[1，2，3]，這些設(shè)備會對時間相關(guān)的加速度載荷產(chǎn)生響應(yīng)，并引起其盛裝的液體介質(zhì)晃蕩，液體隨即又以動壓力的形式反作用于這些設(shè)備，從而改變設(shè)備的振動和變形狀態(tài)。液體晃蕩的作用不僅在設(shè)備自身振動時產(chǎn)生，且振動消失后仍將持續(xù)一段時間。設(shè)備振動時導(dǎo)致的液體晃蕩，不僅受容器本身的振動控制，還會受到慣性和重力的作用[2]；液體晃蕩會對限制其運動的固體表面通過動壓力的方式產(chǎn)生晃蕩力，從而反過來影響容器的運動。顯然，與不含液體介質(zhì)的特種設(shè)備相比，此類設(shè)備在加速度載荷下的結(jié)構(gòu)運動特性會發(fā)生一定的改變，與此同時受到結(jié)構(gòu)振動和液體晃蕩的影響，含液體結(jié)構(gòu)在地震或其他載荷激勵作用下的動力響應(yīng)會非常復(fù)雜。

受加速度載荷激勵引起的液體晃蕩附加力對于設(shè)備安全有很大的影響。當?shù)卣鸢l(fā)生時，如儲罐、球罐、水塔等，地震作用導(dǎo)致的水體晃動而產(chǎn)生的晃蕩力會是導(dǎo)致其結(jié)構(gòu)破壞的主要原因之一[1]，其中破壞形式主要為支座破壞，從而導(dǎo)致設(shè)備傾倒；但同時液體晃蕩可以吸收振動能量，達到減震的目的，如高層建筑的減震水箱[2]，為了滿足使用功能要求，其容量高達上百噸，可能會增加結(jié)構(gòu)的慣性力導(dǎo)致結(jié)構(gòu)震害增加[3]，但若設(shè)計合理水箱可作為調(diào)頻液體阻尼器（Tuned Liquid Damper, TLD）的一部分，利用液體晃蕩過程中產(chǎn)生的動側(cè)壓力提供減振力[4]。因此，從保障特種設(shè)備安全及震后設(shè)備安全評價角度看，為了控制或利用液體晃蕩，優(yōu)化特種設(shè)備抗震結(jié)構(gòu)，充分利用水體晃蕩力來削減地震產(chǎn)生的慣性力，抑制地震對設(shè)備產(chǎn)生破壞，十分有必要對振動作用下液體的動力響應(yīng)進行研究。

液體在振動作用下的動力響應(yīng)由液體晃蕩和液體晃蕩力與支撐結(jié)構(gòu)動力耦合作用兩部分組成，而如何精確獲取液體的晃蕩及其產(chǎn)生的晃蕩力是其中的重點和難點。液體晃蕩的研究方法大致可以歸為理論研究、實驗研究和數(shù)值計算三大類[5-9]。如Chen等[5]通過直接數(shù)值求解歐拉方程的理論方法求解儲罐大幅晃蕩過程；祁江濤[6]和劉文夫等[7]通過VOF方法結(jié)合動網(wǎng)格理論對LNG船液艙晃蕩進行了研究；大連理工大學(xué)岳前進課題組[8]通過搭建大型振動臺對LNG儲液艙的晃蕩載荷分布規(guī)律進行了研究。目前理論研究對復(fù)雜構(gòu)體或粘性較大流體的晃蕩難以獲得精確的解析解；通過振動臺實驗可以模擬地震載荷下儲罐水體晃蕩，但投資巨大，且存在尺度效應(yīng)；隨著計算機技術(shù)的高速發(fā)展，數(shù)值計算方法越來越受到重視，已經(jīng)逐漸成為現(xiàn)代研究液體晃蕩問題的重要手段[10-13]。

VOF（Volume of Fluid）方法[10，11，13]是一種成熟的自由液面追蹤的歐拉數(shù)值方法，它以追蹤流體體積所占網(wǎng)格單元體積分率的途徑描述自由表面變化。本文基于VOF方法對中國特種設(shè)備檢測研究院載荷響應(yīng)實驗室充液球罐晃蕩實驗進行數(shù)值仿真，仿真所得球罐晃蕩附加力與實驗吻合良好，驗證了VOF方法在特種設(shè)備晃蕩研究中的可行性，為后續(xù)研究打下基礎(chǔ)。

1 VOF方法

本文利用ANSYS Fluent?17.0軟件，采用VOF方法模擬水體在球罐中的晃蕩過程，并結(jié)合分段線性化界面計算方法（Piecewise-liner Interface Calculation,PLIC）確定流體間的界面并捕捉氣泡的存在。湍流模型使用了k-ε模型封閉湍流粘度項。在振動過程中水和空氣的總體積基本保持不變，密度幾乎不變，流體的壓縮性對晃蕩作用的影響可忽略，故采用不可壓縮模型[2]；認為晃蕩過程為等溫過程不涉及能量輸運方程的計算；本模擬的兩相流動中的介質(zhì)為空氣和水，不存在反應(yīng)，因此不存在質(zhì)量源相和相間交換項。則根據(jù)以上模型簡化，建立主要的控制方程如下：

體積分率輸運方程：

體積分率之和為1：

式中：

f ——體積分率；

ρ ——密度；

v ——速度；

q ——第q相。

求解過程中僅求解n-1相的輸運方程，而主相的體積分率由體積分率的限制條件即體積分率之和為1確定。

VOF方法中僅求解一個動量輸運方程，認為所有相的速度都是相同的，這種假設(shè)的引入在相界面處速度差異較大的情況下會造成較大的誤差。

動量輸運方程：

式中ρ，μ為控制體內(nèi)的流體密度和粘度，其表達式為：

式中：

μt——湍流粘度，本研究中選用k-ε模型對此項進行封閉；

g ——重力加速度；

F ——外部體積力。

2 球罐振動實驗與模擬

楊智榮等[14]利用中國特種設(shè)備檢測研究院載荷響應(yīng)實驗室的多通道動態(tài)載荷實驗臺對不同充裝液位球罐在不同頻率、振幅外載激勵下的響應(yīng)規(guī)律進行了研究。

實驗系統(tǒng)如圖1、圖2所示，實驗用球罐直徑為1m，液體為水并添加染色劑便于數(shù)字攝影分析。本次實驗的振動輸出為正弦曲線輸出，振幅15mm，頻率1.5Hz；通過加速度傳感器和三維力傳感器測量振動過程中的加速度和球罐水平受力，傳感器分布見圖2。其中傳感器1和3測量主振動方向加速度（X方向），傳感器2測量垂直主振動方向加速度（Y方向）。傳感器4、傳感器5、傳感器6、傳感器7分別安裝在四個支腿底部，可測量四個支腿在X和Y方向的受力，通過加和支腿上的受力再去除裝置的慣性力可以得到晃蕩附加力。

圖1 球罐振動實驗裝置與球罐數(shù)值模擬網(wǎng)格劃分

圖2 球罐水體晃動實驗測點分布

本文選取三個不同液位的實驗過程進行了模擬，盛裝量分別為30%、50%和70%，采用VOF方法耦合k-ε湍流模型，通過用戶自定義函數(shù)（User De fine Function，UDF）添加動量源相的方式實現(xiàn)振動的模擬，輸入的加速度為傳感器1和2的實測加速度，其中測得的加速度（50%液位）如圖3所示。CFD網(wǎng)格采用ANSYS ICEM軟件劃分，全部為六面體網(wǎng)格，總共約245000網(wǎng)格單元，經(jīng)網(wǎng)格無關(guān)性測試表明此網(wǎng)格分辨率可以滿足模擬需求。模擬采用非穩(wěn)態(tài)模擬，時間步長為0.001s，計算時長12s，振動持續(xù)10.5s。罐體表面設(shè)置為無滑移壁面邊界條件。

圖3 裝載量為50%液位時測得的X和Y方向的加速度

表1 球罐晃蕩數(shù)值模擬的物性參數(shù)和模型設(shè)置總結(jié)

ANSYS Fluent中詳細的物性和模型等的設(shè)置見表1。求解在中國特種設(shè)備檢測研究院數(shù)值仿真實驗室的高性能圖形工作站上進行，采用20線程并行計算，單次模擬耗時6h左右。

3 模擬結(jié)果與討論

圖4展現(xiàn)了50%液位時球罐晃蕩模擬不同時刻的液面演變，可以看到在晃蕩初期液體晃蕩速度較小，液面整體傾斜保持基本平直；隨著振動時間增加，速度增大，液面會呈現(xiàn)出波浪狀推進；且在波浪的推進過程中會出現(xiàn)液體與液面的分離拍打，形成氣泡等現(xiàn)象。

圖4 50%液位球罐晃蕩模擬不同時刻的液面演變

不同液位沿振動方向（X方向）的液體晃蕩所引起的附加力與實驗對比，如圖5所示，F(xiàn)orce_X_exp和Force_X_cal分別代表實驗和模擬得到的晃蕩附加力。從圖5中可以看到隨著液位的升高，液體晃蕩所引起的附加力逐漸增大，但并非隨質(zhì)量線性增加；當液面從50%漲到70%時相比液面從30%漲到50%，雖然球罐內(nèi)增加的液體質(zhì)量一致，但晃蕩附加力的增加幅度明顯加大；結(jié)果顯示附加力大小和頻率與實驗值吻合較好，證明此模型可用于振動過程中的球罐結(jié)構(gòu)分析和內(nèi)部晃蕩頻率預(yù)測；在實驗中發(fā)現(xiàn)當液位較低時，水體晃蕩過程中湍動強度較大，產(chǎn)生了較多的氣泡和水花[14]，但在模擬中由于湍流模型的缺陷和VOF方法均勻化的限制，并未很好地捕捉到這一現(xiàn)象，這可能是造成水位較低（30%液面）時實驗與模擬結(jié)果雖然仍在可接受范圍內(nèi)(＜15%)，但明顯大于水位較高（70%液面）時的誤差。在未來的研究中將更多地去考察湍流帶來的影響，以及如何在模擬中準確的進行描述。振動輸入在10.5s處停止，在圖中此時刻之后僅顯示模擬結(jié)果，可以看到在X方向上晃蕩附加力在振動源相消失之后，短時間內(nèi)仍會保持相似的幅值與頻率。

圖5 不同液位沿振動方向（X方向）液體晃蕩模擬得到的附加力與實驗結(jié)果對比

圖6 不同液位沿垂直振動方向（Y方向）液體晃蕩模擬得到的附加力與實驗結(jié)果對比

不同液位沿垂直振動方向（Y方向）的液體晃蕩所引起的附加力與實驗對比，如圖6所示，F(xiàn)orce_Y_exp和Force_Y_cal分別代表實驗和模擬得到的晃蕩附加力。與X方向受力相比，由于振動源的高頻特性，其模擬結(jié)果也呈現(xiàn)高頻狀態(tài)，實驗和模擬結(jié)果基本吻合，實驗測得的液體晃蕩力略大于模擬值；Y方向晃蕩附加力的幅值與X方向相比顯然是不容忽視的，由此可見在實際工程問題中不可忽略垂直主振動方向的液體晃蕩附加力。與X方向類似，隨著液位升高，液體晃蕩所引起的附加力逐漸增大，但并非隨質(zhì)量線性增加，從50%到70%液面的晃蕩附加力增幅明顯更大。在振動源相消失后，Y方向的晃蕩附加力表現(xiàn)出與X方向附加力不同的響應(yīng)規(guī)律，其幅值會立即下降，這可能是由于Y方向的液體晃蕩呈現(xiàn)高頻的特點造成的。圖5中也能體現(xiàn)這一規(guī)律，圖3中X方向的振動輸入加速度，可以認為是由兩個波疊加而生成的，一個幅值大頻率低的長波和一個類似于Y方向加速度的高頻短波，圖5中停止振動后液體晃蕩力的曲線立即變得平滑，可以認為是高頻短波的影響在振動消失之后很快衰減。

為考察液體晃蕩附加力的頻率是否與振動頻率一致，對模擬得到的晃蕩附加力進行了正弦函數(shù)擬合，例如50%液位的擬合結(jié)果見圖7。振動源的頻率為1.5Hz，而30%、50%和70%液面通過擬合分析得到的頻率分別為1.497Hz、1.485Hz和1.503Hz，可見液體晃蕩產(chǎn)生的附加力的頻率與振動頻率基本一致。

圖7 50%液面時模擬得到的晃蕩附加力正弦函數(shù)數(shù)值擬合結(jié)果

4 結(jié)論與展望

本文通過VOF方法耦合界面重構(gòu)模型和k-ε湍流模型對不同液位下的球罐晃蕩行了模擬，對球罐壁面壓力進行積分得到了晃蕩引起的附加力。模擬過程實現(xiàn)了實際測量加速度與計算流體軟件的耦合，有利于本模型的推廣。根據(jù)主振動方向和垂直其方向的晃蕩附加力的實驗和模擬結(jié)果，發(fā)現(xiàn)在考察球罐抗震安全性時垂直振動方向的高頻振動的影響不可忽略；但液面較低時（30%液面）湍流強度較大，模擬與實驗結(jié)果的偏差較大，湍流模型在球罐晃蕩模擬中的作用需要進一步研究；從對比結(jié)果看出VOF方法可以較好預(yù)測振動狀態(tài)下的液體晃蕩附加力，證明該模型可用于大型球罐振動晃蕩過程的研究。未來可通過模擬對大型儲液設(shè)備的晃蕩附加力的影響因素進行更加深入的研究，并引入流固耦合仿真對設(shè)備進行含液體晃蕩力的結(jié)構(gòu)分析，為儲液設(shè)備達到本質(zhì)安全發(fā)揮重要的作用。

參考文獻

[1]Sezen H， Livaoglu R， Dogangun A．Dynamic Analysis and Seismic Performance Evaluation of Above-ground Liquid-containing Tanks [J]．Engineering Structure， 2008，30(03)： 794-803．

[2]尚春雨．帶自由水面的結(jié)構(gòu)地震反應(yīng)分析方法及應(yīng)用[D]．上海：上海交通大學(xué)，2008．

[3]劉云賀，張俊發(fā)，田潔，等．頂部設(shè)置水箱的建筑結(jié)構(gòu)控震研究[J]．西安理工大學(xué)學(xué)報，1999，15(04)：71-75．

[4]榮嶠．高層結(jié)構(gòu)利用調(diào)諧液體阻尼器（TLD）控制水平地震反應(yīng)的研究[D]．北京：北京交通大學(xué)，2009．

[5]Chen W， Haroun M A， Liu F．Large amplitude liquid sloshing in seismically excited tanks [J]．Earthquake Engineering and Structural Dynamics， 1996， 25： 653-669．

[6]祁江濤，顧民．LNG船液艙晃蕩的數(shù)值模擬[J]．中國造船，2007（b11)： 541-549．

[7]劉文夫，薛鴻祥，唐文勇，等．獨立B型LNG棱形液艙晃蕩載荷數(shù)值分析[J]．船舶工程，2015，37(07)：22-25+72．

[8]衛(wèi)志軍，阮詩倫，岳前進，等．超大型儲液艙內(nèi)晃蕩載荷實驗研究[J]．海洋工程裝備與技術(shù)，2014，1(01)：55-61．

[9]衛(wèi)志軍，岳前進，張文首，等．大尺度儲艙液體晃蕩砰擊壓力測量方法研究[J]．中國科學(xué)：物理學(xué)力學(xué) 天文學(xué)，2014，44（07）：746-758．

[10]ANSYS Inc．， ANSYS FLUENT 17．0 Documentation [Z]．2016．

[11]Youngs D L．Time-Dependent Multi-Material Flow with Large Fluid Distortion， Numerical Methods for Fluid Dynamics [C]．Morton K W and Baines M J， editors．1982， Academic Press．

[12]Hasheminejad M S， Mohammadi M M， Jarrahi M．Liquid Sloshing in Partly-filled Laterally-exited Circular Tanks Equipped with Baffles [J]．Journal of Fluids and Structures， 2014， 44： 97-114．

[13]龔國毅，趙成璧，唐友宏，等．基于VOF和浸入邊界法的黏性二相流模型對LNG液艙晃蕩的數(shù)值模擬[J]．科學(xué)技術(shù)與工程，2013，13(10)：2630-2634．

[14]楊智榮．大型球罐結(jié)構(gòu)抗震安全保障中的關(guān)鍵技術(shù)研究[R]．北京：2016．