壁面結(jié)構(gòu)形狀對流體晃蕩影響的數(shù)值研究

衛(wèi)志軍，杜祥璞，鄒國良，吳錘結(jié)，翟鋼軍

（1.大連理工大學(xué)a.海岸和近海工程國家重點實驗室，b.航空航天學(xué)院，遼寧 大連 116023；2.南京水利科學(xué)研究院水文水資源與水利工程科學(xué)國家重點實驗室，南京 210029）

0 引 言

中國的海洋油氣資源儲量豐富，向南海石油開發(fā)進(jìn)軍，是我國對南海海域?qū)嵭泄苤频膽?zhàn)略措施，還可以增加我國的石油產(chǎn)量和巨大的經(jīng)濟(jì)利益。浮式液化天然氣（Floating Liquified Natural Gas，F(xiàn)LNG）終端具有開采周期短、可獨(dú)立作業(yè)、可根據(jù)實際開采油氣田生產(chǎn)狀況靈活配置、無需管道輸送、可回收和可運(yùn)移等特點，非常適合我國南海深遠(yuǎn)海油氣田開發(fā)[1]。特殊浮式開采平臺的噸位和液艙容積越來越大，且無任何載液率限制，這不僅影響船舶的運(yùn)動[2-3]，導(dǎo)致流體對儲艙結(jié)構(gòu)產(chǎn)生較大的晃蕩甚至砰擊載荷，進(jìn)一步引起艙壁結(jié)構(gòu)強(qiáng)度失效，甚至浮式平臺失穩(wěn)[4-8]。因此如何合理有效地抑制儲艙內(nèi)流體晃蕩是亟待解決的工程科學(xué)問題。

國內(nèi)外學(xué)者針對如何抑制流體晃蕩開展了一系列研究，主要采用在液艙內(nèi)增設(shè)阻晃隔板或者在自由液面加設(shè)浮式裝置的方式。在液艙內(nèi)加阻晃隔板方面，Akyildiz等[9]研究了不同位置的壓力分布以及3D效應(yīng)對液體晃動的影響。薛米安等[10]通過實驗在液艙中間加不同形式隔板并監(jiān)測艙壁的沖擊壓力。Kandasamy等[11]數(shù)值模擬研究了三維液艙不同填裝水平下，橫向、斜向、常規(guī)、局部縱向和橫向加速度作用下?lián)醢逶O(shè)計的有效性。衛(wèi)志軍等[12]通過模型試驗研究了近二維淺水條件下，液艙在大振幅簡諧側(cè)向激振作用和高固體率的垂直阻晃隔板作用下的砰擊現(xiàn)象。Singal等[13]用數(shù)值模擬方法對煤油液體界面進(jìn)行瞬態(tài)模擬，發(fā)現(xiàn)在燃料箱中引入擋板后，燃料箱內(nèi)的晃動明顯減小。對有內(nèi)結(jié)構(gòu)和無內(nèi)結(jié)構(gòu)的儲罐進(jìn)行了數(shù)值研究。蔡忠華等[14]數(shù)值討論了強(qiáng)肋框、水平桁及制蕩艙壁對于壓力及液艙內(nèi)液體固有頻率的影響。在自由液面加設(shè)浮式裝置方面，寧德志等[15]通過實驗對三維液艙內(nèi)柱型浮子式減晃結(jié)構(gòu)的水動力特性開展了研究。Koh等[16]利用改進(jìn)后的CPM成功地模擬了帶約束浮動擋板（CFB）的水晃動。Yu等[17]在六自由度運(yùn)動平臺上的剛性矩形液艙內(nèi)進(jìn)行了不同固體率的雙穿孔浮板抗晃動模型的掃頻試驗。Jithu等[18]基于RANSE數(shù)值模擬研究了有球型浮動擋板和無球型浮動擋板情況下的儲罐結(jié)構(gòu)變形和儲罐壁面壓力。李雪艷等[19]考慮波面高程、透射系數(shù)、反射系數(shù)、能耗系數(shù)及耗散波能等指標(biāo)，對平板式、上弧板式和下弧板式三種透空堤開展了消浪性能數(shù)值研究。

由于上述兩種抑制方式均對工藝要求甚高，在實際建造過程中較難實現(xiàn)[15,20-21]。真實船載LNG液艙的建造中，法國GTT公司設(shè)計的液化天然氣儲艙MARKIII表面具有凸起結(jié)構(gòu)。從壁面結(jié)構(gòu)自身結(jié)構(gòu)特性開展流體晃蕩抑制可以滿足建造工藝，有利于工程實現(xiàn)。然而國內(nèi)外針對壁面凸起結(jié)構(gòu)抑制流體晃蕩影響的研究較少。Graczyk等[22]通過實驗方法分析在二維模型低液位的工況下凸起結(jié)構(gòu)對流體晃蕩荷載的影響，研究表明凸起物使流體沖擊壁面時引入大量空氣，并且水平凸起物下面的壓力比上面的更大。Moirod等[23]采用數(shù)值方法對比分析一個光滑的楔形體和一個凸起結(jié)構(gòu)的楔形體的入水砰擊試驗，結(jié)果表明凸起結(jié)構(gòu)可以增大晃動壓力，建議設(shè)計超大型船載儲艙時，需考慮凸起結(jié)構(gòu)對砰擊的影響。目前尚未有針對壁面凸起結(jié)構(gòu)抑制流體晃蕩的系統(tǒng)研究。為此，本文采用數(shù)值方法，從物理機(jī)理角度系統(tǒng)分析三種不同形狀凸起結(jié)構(gòu)對流體晃蕩的抑制效果。結(jié)果表明，合理設(shè)計凸起結(jié)構(gòu)能夠改變流體砰擊模式，有效降低晃蕩荷載沖量，進(jìn)而有效抑制流體晃蕩。該研究結(jié)果能為儲艙抗砰擊設(shè)計提供科學(xué)參考。

1 計算模型參數(shù)

矩形液艙模型尺寸及波高與壓力測點如圖1所示，模型尺寸長（L）1 m、高（H）1 m、寬（B）0.1 m，壁面凸起塊凸起高度0.05 m，與艙壁接觸面為邊長0.05 m的正方形，采用ICEM軟件建模、結(jié)構(gòu)化網(wǎng)格和FLUENT計算。通過文獻(xiàn)調(diào)研得知，在流體晃蕩過程中氣體的氣墊效應(yīng)是不可忽略的因素[24-25]。因此本文考慮氣體可壓縮性，采用VOF模型對自由液面進(jìn)行捕捉[26-28]，湍流模型選擇Standardk-ε模型，計算參數(shù)設(shè)置如表1所示。

表1 數(shù)值模擬參數(shù)設(shè)置Tab.1 Numerical simulation parameter settings

圖1 矩形液艙尺寸示意圖（單位：mm）Fig.1 Schematic diagram of rectangular tank size(Unit:mm)

2 網(wǎng)格無關(guān)性檢驗

本文采用二維模型，分別針對5種網(wǎng)格進(jìn)行無關(guān)性檢驗，參數(shù)設(shè)置如表2所示，時間步長為0.001 s。考慮到淺水工況下，流體運(yùn)動較劇烈，因此設(shè)計h/H=0.125，h為載液高度。激勵形式采用單自由度正弦激勵，

表2 網(wǎng)格無關(guān)性計算工況Tab.2 Grid independence calculation condition

式中，A為激勵幅值，f為激勵頻率，t為激勵時間。網(wǎng)格無關(guān)性檢驗中，A為0.1m，f/f0=0.74，f0為自由液面最低階固有頻率。通過對比晃蕩荷載和近壁面波高發(fā)現(xiàn)，當(dāng)網(wǎng)格增加到工況3時，計算結(jié)果已經(jīng)趨于收斂，如圖2所示。

圖2 網(wǎng)格無關(guān)性檢驗Fig.2 Grid independence test

進(jìn)一步分析計算誤差，

式中，η為壓力峰值或波高峰值，i=1，2，3，4。分析結(jié)果見表3。

表3 誤差分析Tab.3 Error analysis

綜合考慮計算時間與精度，二維情況下工況3網(wǎng)格間距滿足計算要求，因此本文選擇網(wǎng)格間距為5 mm開展后續(xù)研究。

3 數(shù)值可靠性驗證

為了保證數(shù)值方法的可靠性，本項目在大連理工大學(xué)工業(yè)裝備及結(jié)構(gòu)分析國家重點實驗室開展了相關(guān)模型試驗研究，將數(shù)值結(jié)果與實驗結(jié)果進(jìn)行系統(tǒng)對比，校核數(shù)值模型。實驗液艙模型見圖3（a），壓力測點見圖1。由于流體非線性影響[29-30]，結(jié)果的峰值稍有差別，但在晃蕩荷載趨勢和量級方面吻合較好，如圖3（b）所示，因此本次數(shù)值方法可用于后續(xù)研究工作。

圖3 數(shù)值模擬與實驗壓力結(jié)果對比，h/H=0.125，f/f0=0.74Fig.3 Comparison of numerical simulation and experimental pressure results,h/H=0.125,f/f0=0.74

4 結(jié)果與討論

本文屬于概念設(shè)計初步研究工作，僅在二維情況下討論單一壁面結(jié)構(gòu)對流體晃蕩的影響。考慮到流體晃蕩的能量主要分布在自由液面最低階固有頻率附近，因此數(shù)值試驗的無量綱激勵頻率f/f0分別為0.8、1和1.2。激勵函數(shù)如式（1），載液率為0.125，激勵振幅為0.1 m。本文以三種不同壁面結(jié)構(gòu)的形狀（三角形、梯形、正方形）為變量，并以光滑壁面作為對照試驗，從時間和空間分布方面系統(tǒng)分析壁面結(jié)構(gòu)形狀對晃蕩荷載的影響，最后分析其對波高的影響。

4.1 壁面結(jié)構(gòu)形狀對晃蕩荷載時間分布的影響

為了細(xì)致分析壁面結(jié)構(gòu)形狀對流場的影響，根據(jù)凸起結(jié)構(gòu)的位置，將左側(cè)艙壁分成上下兩個區(qū)域，如圖4所示，分別從上和下區(qū)域各選擇一個具有代表性的測點P2和P3。在最低階固有頻率附近，具有四種不同壁面結(jié)構(gòu)形式的液艙中壓力時程曲線對比如圖5所示。

圖4 上部區(qū)域和下部區(qū)域Fig.4 Schematic diagram of upper and lower regions

圖5 壓力時程對比Fig.5 Comparision of the pressure time histories

對于下部區(qū)域P2測點，流體沖擊具有壁面凸起結(jié)構(gòu)時，壁面晃蕩壓力表現(xiàn)為雙峰，如圖6所示，且流體砰擊的第一個峰值較尖銳，持續(xù)時間較短。此外，不同形狀的凸起結(jié)構(gòu)下的壓力均大于光滑壁面壓力：當(dāng)激勵頻率f/f0=0.8時，梯形凸起結(jié)構(gòu)壓力峰值最大，約為相同工況下光滑壁面壓力的5倍；當(dāng)激勵頻率f/f0=1時，三角形凸起結(jié)構(gòu)壓力峰值最大，約為光滑壁面的1.5倍；當(dāng)激勵頻率f/f0=1.2時，梯形凸起結(jié)構(gòu)壓力峰值最大，約為光滑壁面3倍。產(chǎn)生上述現(xiàn)象的原因?qū)⒃诤罄m(xù)給出。然而，光滑壁面工況中由于重力引起的第二個壓力峰值均大于具有壁面結(jié)構(gòu)工況中的壓力，且隨著頻率的增大，差距越顯著。分析認(rèn)為由于凸起結(jié)構(gòu)的存在，使得流體能量損耗增加，流體動能轉(zhuǎn)化為勢能減少，流體沿著壁面爬高降低，進(jìn)而導(dǎo)致流體晃蕩的第二個壓力峰值下降。此外，在流體沿著壁面向上爬升過程中，需繞過壁面凸起結(jié)構(gòu)作用到上部區(qū)域，這一過程導(dǎo)致同一周期內(nèi)有壁面凸起結(jié)構(gòu)時的壓力峰值比光滑壁面在時間上發(fā)生得晚，導(dǎo)致壓力時程曲線具有一定的相位差。

圖6 壓力雙峰曲線Fig.6 Pressure bimodal curve

從上部區(qū)域P3測點的峰值看，受凸起結(jié)構(gòu)的影響，流體對上部區(qū)域作用力均有所下降，在f/f0=1和f/f0=1.2時降低較為明顯。所以對于下部區(qū)域，壁面凸起結(jié)構(gòu)引起測點壓力第一個峰值上升，第二個峰值下降；對于上部區(qū)域，壁面結(jié)構(gòu)均會促使流體對艙壁的作用力顯著下降。

通過壓力時程對比，為了深入分析壁面結(jié)構(gòu)抑制流體晃蕩的物理機(jī)理，針對半個周期內(nèi)自由液面和晃蕩荷載的演化過程，針對代表性的具有矩形凸起結(jié)構(gòu)的壁面與光滑壁面工況開展細(xì)致分析。當(dāng)f/f0＝0.8時，光滑壁面液艙半個周期內(nèi)的流場隨時間變化如圖7所示，液艙從原點向左運(yùn)動，液艙內(nèi)一列行波向左側(cè)壁面?zhèn)鞑ィ▓D7（a）#1），該行波具有顯著的水頭，隨著左側(cè)壁面附近波谷向上抬升，水頭高度逐漸減低，流體直接沖擊到壁面約0.08 m處（圖7（a）#2），凸起下部區(qū)域形成了由于砰擊產(chǎn)生的第一個峰值（圖7（c）#2），且此時流體尚未到達(dá)上部區(qū)域（圖7（d）#2）。隨著液艙運(yùn)動到左側(cè)最大位移處，液艙運(yùn)動速度逐漸降低到0，由于慣性力的作用，流體沿壁面向上爬升（圖7（a）#3）直至動能全部轉(zhuǎn)化為勢能，流體運(yùn)動到垂向最大位置。隨后，流體在重力作用下，開始回落，即勢能逐漸轉(zhuǎn)化為動能，下部區(qū)域形成壓力雙峰曲線的第二個峰值（圖7（c）#4），而由于大部分流體已離開上部區(qū)域，因此上部區(qū)域在重力作用下的流體作用力較?。▓D7（d）#4）。當(dāng)液艙反向運(yùn)動，流體也向右側(cè)艙壁傳播，如圖7（a）#5和#6所示。

圖7 半周期內(nèi)，f/f0＝0.8工況下光滑壁面流場變化Fig.7 Flow evalution during half cycle for smooth wall surface subjected to f/f0=0.8

考慮到對流體晃蕩的抑制效果，本文選擇效果較好的正方形凸起壁面流場隨時間變化進(jìn)行分析，液艙位移與半個周期內(nèi)的流場隨時間變化如圖8所示，液艙從原點向左運(yùn)動，液艙內(nèi)一列行波向左側(cè)壁面?zhèn)鞑ィ▓D8（a）#1）。當(dāng)流體沖擊到左側(cè)壁面被凸起結(jié)構(gòu)分割成兩部分時，下層流體受阻擋被困在液艙底面與凸起結(jié)構(gòu)之間，并與凸起結(jié)構(gòu)圍成一個大的封閉氣腔（圖8（a）#2），流體通過壓縮空氣作用在壁面上，形成壓力雙峰曲線第一個峰值（圖8（c）#2），正是由于形成氣腔導(dǎo)致凸起結(jié)構(gòu)第一個壓力峰值比光滑壁面第一個壓力峰值提前約0.02 s，這說明壁面凸起結(jié)構(gòu)可以延長下部區(qū)域荷載作用時間。隨著液艙速度下降，上層流體沖擊到正方形壁面凸起結(jié)構(gòu)，其運(yùn)動方向和自由液面形狀發(fā)生了比較大的變化。水頭躍起到半空（圖8（a）#3）（三角形和梯形凸起結(jié)構(gòu)下流體產(chǎn)生翻卷現(xiàn)象），此時壓力在雙峰之間波谷處（圖8（c）#3）以靜水壓力為主，而發(fā)生水躍將直接導(dǎo)致上部區(qū)域壁面凸起結(jié)構(gòu)壓力峰值比光滑壁面在時間上發(fā)生得晚。在（圖8（a）#4）中液艙速度為0時，流體受重力影響開始回落形成壓力雙峰曲線的第二個峰值（圖8（c）#4）。

圖8 半周期內(nèi)，f/f0＝0.8工況下正方形凸起壁面流場變化Fig.8 Flow evalution during half cycle for square raised wall surface subjected to f/f0=0.8

通過對流場沖擊過程中自由液面物理現(xiàn)象的分析，發(fā)現(xiàn)添加壁面結(jié)構(gòu)改變產(chǎn)生最大砰擊荷載的原因：流體直接砰擊光滑側(cè)壁，產(chǎn)生最大的砰擊荷載（圖7（a）#2），轉(zhuǎn)為氣體，在流體擠壓作用下，轉(zhuǎn)而作用具有壁面結(jié)構(gòu)的側(cè)壁，產(chǎn)生最大砰擊荷載（圖8（a）#2）。針對圖7#2和圖8#2流場瞬態(tài)空氣密度分布（圖9）開展分析，發(fā)現(xiàn)空氣密度分布在水平方向上越接近左側(cè)壁面，密度越大；豎直方向上越接近自由液面密度越大，在凸起結(jié)構(gòu)下側(cè)與壁面接觸點密度達(dá)到最大。

圖9 不同形狀凸起結(jié)構(gòu)空氣密度分布Fig.9 Air density distribution of convex structure with different shapes

當(dāng)流體到達(dá)壁面，受慣性作用和豎直液艙壁的影響有向上爬升的趨勢，光滑壁面沖擊位置附近空氣擴(kuò)散和流體向上爬升不受阻礙。相反地，由于凸起結(jié)構(gòu)阻礙了流體向上爬升，即凸起邊界對流體產(chǎn)生了擠壓作用。由于空氣可壓縮性，從空氣密度云圖來看，相對于光滑壁面，壁面結(jié)構(gòu)沖擊位置處氣體密度較大，說明空氣被壓縮，同理液體也受擠壓。此時流體受擠壓的力會反作用于壁面，從而導(dǎo)致該區(qū)域的壁面壓力增加。從另一個角度分析，空氣被局部壓縮必然導(dǎo)致該區(qū)域氣體壓力增加，相當(dāng)于凸起結(jié)構(gòu)下部區(qū)域受到了較大的局部氣壓，該局部區(qū)域的絕對氣壓大于整個液艙所處環(huán)境的氣壓。因此，無論是空氣和水直接受到向上的阻礙還是空氣被壓縮，都是壁面凸起結(jié)構(gòu)的下壁面壓力增加的原因。

當(dāng)砰擊發(fā)生后，壁面凸起結(jié)構(gòu)對流場湍動能的影響如圖10所示。湍動能主要集中于凸起結(jié)構(gòu)附近，即凸起結(jié)構(gòu)是導(dǎo)致流場湍流增強(qiáng)的主要原因。湍動能的增加引起了動能耗散的增加如圖11所示，添加壁面凸起結(jié)構(gòu)時湍動能耗散率顯著提高，所以砰擊發(fā)生后，由于結(jié)構(gòu)的存在導(dǎo)致能量消耗，從而使得流體晃蕩荷載降低。這是導(dǎo)致凸起結(jié)構(gòu)下部區(qū)域第二峰值和上部區(qū)域壓力相對于光滑壁面下降的主要原因。

圖10 不同形狀凸起結(jié)構(gòu)湍動能Fig.10 Turbulence kinetic energy of convex structure with different shapes

圖11 不同形狀凸起結(jié)構(gòu)湍動能耗散率Fig.11 Turbulent dissipation rate of convex structure with different shapes

4.2 壁面結(jié)構(gòu)形狀對晃蕩荷載空間分布的影響

4.2.1 晃蕩壓力空間分布

考慮到晃蕩荷載的不確定性，本文采用5個周期壓力峰值的平均值分析壁面結(jié)構(gòu)形狀對晃蕩壓力空間分布的影響，如圖12所示。隨著激勵頻率的增加，流體沖擊壁面的位置逐漸提高。隨著沖擊位置的變化，凸起結(jié)構(gòu)的形狀對荷載的影響敏感度較大。從整體來看凸起結(jié)構(gòu)對上下兩個區(qū)域影響結(jié)果相反，即凸起結(jié)構(gòu)會引起下部區(qū)域產(chǎn)生較大的局部壓力，而上部區(qū)域的壓力則相對于光滑壁面有所下降。

圖12 晃蕩壓力峰值平均值空間分布Fig.12 Spatial distribution of the mean value of the sloshing pressure peak

4.2.2 晃蕩荷載沖量空間分布

考慮到液艙的安全性，結(jié)構(gòu)響應(yīng)是不可忽略的，相比于瞬間較大的局部壓力，持續(xù)的壓力作用，即晃蕩荷載沖量，才是影響結(jié)構(gòu)響應(yīng)的重要因素。本文通過計算沖量I評價不同液艙中砰擊壓力的動力特性和其空間分布，

式中，t1和t2分別代表砰擊前、后的時間，p是由于砰擊產(chǎn)生的局部壓力。

沖量為壓力曲線下的紅色陰影面積，如圖13（a）所示。當(dāng)f/f0=0.8時，壁面凸起結(jié)構(gòu)對下部區(qū)域影響不明顯，導(dǎo)致上部區(qū)域沖量下降較大；當(dāng)f/f0=1時，壁面沖量相對光滑壁面下降均較明顯，在上部區(qū)域沖量甚至下降一半以上；當(dāng)f/f0=1.2時，壁面凸起結(jié)構(gòu)工況中下部區(qū)域沖量下降明顯，上部區(qū)域梯形和正方形凸起結(jié)構(gòu)可以引起沖量下降，而三角形凸起結(jié)構(gòu)超過了光滑壁面沖量。從這三種不同頻率工況來看，激勵頻率f/f0=1時，壁面結(jié)構(gòu)可以更有效抑制流體晃蕩。綜上所述，采用壁面結(jié)構(gòu)形狀抑制流體晃蕩是行之有效的方法。但需設(shè)計合理的壁面凸起結(jié)構(gòu)形狀，以有效降低荷載沖量，從而抑制晃蕩。

圖13 沖量空間分布對比Fig.13 Comparison of impulse spatial distributions

4.3 壁面結(jié)構(gòu)對波高的影響

本文選取穩(wěn)定狀態(tài)波高時程數(shù)據(jù)進(jìn)行分析，各個周期波峰值標(biāo)準(zhǔn)差數(shù)量級很小，如表4所示，說明峰值波動較小，可以對峰值進(jìn)行平均分析。差值計算公式為

表4 不同壁面結(jié)構(gòu)的波高平均峰值Tab.4 Average peak of wave height of different wall structures

式中，h0為光滑壁面波高平均峰值，hi為有凸起結(jié)構(gòu)壁面波高平均峰值，i=1,2,3分別代表三角形、梯形和正方形凸起壁面。

不同壁面結(jié)構(gòu)波高曲線趨勢基本一致，波高峰值的第二個峰值大于第一個峰值，如圖14所示。

圖14 自由液面高度對比Fig.14 Comparison of wave heights

本文以一個周期內(nèi)波高最大值，即第二峰值為分析對象。激勵頻率f/f0＝0.8時，壁面凸起結(jié)構(gòu)對波高的抑制比較明顯；激勵頻率f/f0＝1時，壁面結(jié)構(gòu)對波高依然具有一定的抑制效果，但是抑制效果較低頻率有所下降；激勵頻率f/f0＝1.2時，由于加壁面凸起結(jié)構(gòu)自由液面破碎較嚴(yán)重，波高反而增大了。結(jié)果表明，在最低階固有頻率及其以下低頻率，壁面凸起結(jié)構(gòu)對波高有抑制作用，隨著頻率的增加抑制作用有下降趨勢。當(dāng)頻率高于最低階固有頻率時，壁面凸起結(jié)構(gòu)可導(dǎo)致波高的增加。

5 結(jié) 論

本文采用數(shù)值方法，系統(tǒng)分析了不同形狀壁面結(jié)構(gòu)形狀對流體晃蕩的影響，主要結(jié)論如下：

（1）壁面凸起結(jié)構(gòu)改變了流體砰擊模式。

（2）流體從直接砰擊壁面變?yōu)榘鼩馀蓍g接砰擊壁面，因此氣泡的產(chǎn)生影響了晃蕩荷載的時空分布及其沖量和波高，并增強(qiáng)流體能量耗散。由此可見，氣體可壓縮性是研究流體晃蕩需考慮的重要因素之一。

（3）合理設(shè)計壁面結(jié)構(gòu)的形狀可以有效地抑制晃蕩，未來將系統(tǒng)研究不同載液率、凸起結(jié)構(gòu)數(shù)量以及自由液面與凸起結(jié)構(gòu)相對位置對流體晃蕩的影響等問題。