基于VOF法的球罐內液體晃蕩數值模擬及載荷計算

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

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

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

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

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

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

1 VOF方法

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

體積分率輸運方程：

體積分率之和為1：

式中：

f ——體積分率；

ρ ——密度；

v ——速度；

q ——第q相。

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

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

動量輸運方程：

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

式中：

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

g ——重力加速度；

F ——外部體積力。

2 球罐振動實驗與模擬

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

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

圖1 球罐振動實驗裝置與球罐數值模擬網格劃分

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

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

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

表1 球罐晃蕩數值模擬的物性參數和模型設置總結

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

3 模擬結果與討論

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

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

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

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

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

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

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

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

4 結論與展望

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

參考文獻

[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]尚春雨．帶自由水面的結構地震反應分析方法及應用[D]．上海：上海交通大學，2008．

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

[4]榮嶠．高層結構利用調諧液體阻尼器（TLD）控制水平地震反應的研究[D]．北京：北京交通大學，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船液艙晃蕩的數值模擬[J]．中國造船，2007（b11)： 541-549．

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

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

[9]衛志軍，岳前進，張文首，等．大尺度儲艙液體晃蕩砰擊壓力測量方法研究[J]．中國科學：物理學力學 天文學，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液艙晃蕩的數值模擬[J]．科學技術與工程，2013，13(10)：2630-2634．

[14]楊智榮．大型球罐結構抗震安全保障中的關鍵技術研究[R]．北京：2016．