壁面結構形狀對流體晃蕩影響的數值研究

衛志軍，杜祥璞，鄒國良，吳錘結，翟鋼軍

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

0 引 言

中國的海洋油氣資源儲量豐富，向南海石油開發進軍，是我國對南海海域實行管制的戰略措施，還可以增加我國的石油產量和巨大的經濟利益。浮式液化天然氣（Floating Liquified Natural Gas，FLNG）終端具有開采周期短、可獨立作業、可根據實際開采油氣田生產狀況靈活配置、無需管道輸送、可回收和可運移等特點，非常適合我國南海深遠海油氣田開發[1]。特殊浮式開采平臺的噸位和液艙容積越來越大，且無任何載液率限制，這不僅影響船舶的運動[2-3]，導致流體對儲艙結構產生較大的晃蕩甚至砰擊載荷，進一步引起艙壁結構強度失效，甚至浮式平臺失穩[4-8]。因此如何合理有效地抑制儲艙內流體晃蕩是亟待解決的工程科學問題。

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

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

1 計算模型參數

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

表1 數值模擬參數設置Tab.1 Numerical simulation parameter settings

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

2 網格無關性檢驗

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

表2 網格無關性計算工況Tab.2 Grid independence calculation condition

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

圖2 網格無關性檢驗Fig.2 Grid independence test

進一步分析計算誤差，

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

表3 誤差分析Tab.3 Error analysis

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

3 數值可靠性驗證

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

圖3 數值模擬與實驗壓力結果對比，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 結果與討論

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

4.1 壁面結構形狀對晃蕩荷載時間分布的影響

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

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

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

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

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

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

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

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

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

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

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

圖9 不同形狀凸起結構空氣密度分布Fig.9 Air density distribution of convex structure with different shapes

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

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

圖10 不同形狀凸起結構湍動能Fig.10 Turbulence kinetic energy of convex structure with different shapes

圖11 不同形狀凸起結構湍動能耗散率Fig.11 Turbulent dissipation rate of convex structure with different shapes

4.2 壁面結構形狀對晃蕩荷載空間分布的影響

4.2.1 晃蕩壓力空間分布

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

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

4.2.2 晃蕩荷載沖量空間分布

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

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

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

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

4.3 壁面結構對波高的影響

本文選取穩定狀態波高時程數據進行分析，各個周期波峰值標準差數量級很小，如表4所示，說明峰值波動較小，可以對峰值進行平均分析。差值計算公式為

表4 不同壁面結構的波高平均峰值Tab.4 Average peak of wave height of different wall structures

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

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

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

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

5 結 論

本文采用數值方法，系統分析了不同形狀壁面結構形狀對流體晃蕩的影響，主要結論如下：

（1）壁面凸起結構改變了流體砰擊模式。

（2）流體從直接砰擊壁面變為包裹氣泡間接砰擊壁面，因此氣泡的產生影響了晃蕩荷載的時空分布及其沖量和波高，并增強流體能量耗散。由此可見，氣體可壓縮性是研究流體晃蕩需考慮的重要因素之一。

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