衛(wèi)志軍，張文首，王安良，董玉山，胡方源，岳前進(jìn)

（大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116024）

0 引 言

隨著超大型浮式液化天然氣（floating liquid natural gas，F(xiàn)LNG）船的研發(fā)，超大型船載儲液艙內(nèi)的液體晃蕩問題成為FLNG 船設(shè)計中的關(guān)鍵問題.儲液艙內(nèi)液體晃蕩是當(dāng)外激激勵頻率與容器內(nèi)部液體自由液面的固有頻率接近時，液體產(chǎn)生的劇烈共振運動.

研究液體晃蕩問題主要有理論、實驗和數(shù)值3種方法.晃蕩的高度隨機(jī)性和非線性，使得理論方法具有一定局限性［1］.實驗方法在晃蕩問題的研究中仍然是最可靠的方法.不過物理實驗設(shè)備的成本非常昂貴，而且實驗的周期都相對很長，所以實驗方法的經(jīng)濟(jì)性較差.隨著計算機(jī)技術(shù)的快速發(fā)展和計算機(jī)容量的迅速擴(kuò)大，計算機(jī)的計算能力得到很大提升且對物理現(xiàn)象的模擬也更加精確，可以采用數(shù)值方法開展長時間大量的仿真研究.因此，數(shù)值方法可以很好地彌補(bǔ)前二者的缺點，成為研究晃蕩問題的主要方法之一.基于網(wǎng)格的數(shù)值模擬方法已經(jīng)在計算力學(xué)領(lǐng)域被廣泛應(yīng)用，例如有限元法、有限體積法和有限差分法等.液體晃蕩具有強(qiáng)非線性自由液面（波浪破碎和翻卷等），基于網(wǎng)格法的數(shù)值方法在模擬自由液面的強(qiáng)非線性時存在相當(dāng)大的困難且對網(wǎng)格質(zhì)量有很強(qiáng)的依賴性［2］.而光滑粒子流體動力學(xué)（smoothed particle hydrodynamics，SPH）方法是一種無網(wǎng)格Lagrangian算法.該方法于1977年分別由Lucy和Gingold等首次同時提出，起初這種方法主要用來解決星體引力天體物理方面所涉及的相關(guān)力學(xué)問題［3－4］.當(dāng)時，這種方法在精度和穩(wěn)定性方面還存在不足，所以它并未得到廣泛的應(yīng)用.后來Monaghan對SPH 方法進(jìn)行了改進(jìn)，首次將其運用到了模擬潰壩、鉆削，以及造波問題［5］.SPH 方法不僅省去了網(wǎng)格劃分的煩瑣步驟，還能很好地模擬自由液面大變形問題.因此，SPH 方法被廣泛應(yīng)用于爆炸、彈道和晃蕩等問題的研究.

影響儲液艙內(nèi)液體晃蕩的因素很多，如載液率、激勵方式、激勵振幅和激勵頻率等［1，6］.很多學(xué)者采用SPH 方法對儲液艙內(nèi)液體晃蕩開展了研究.陳正云等采用SPH 方法研究二維液艙大振幅非諧振頻率激勵下液體晃蕩的自由液面的波形變化［7］.崔巖等采用SPH 方法研究二維矩形液艙在共振頻率附近縱蕩、縱搖，以及二者綜合效應(yīng)下波高的變化［8］.李大鳴等采用SPH 方法研究不同橫蕩振幅對二維矩形液艙非共振頻率激勵下自由液面的影響［9］.此外，李大鳴等還采用SPH 方法研究遠(yuǎn)離共振頻率的激勵頻率下，液艙內(nèi)中高載液率對自由液面的影響［10］.

上述研究工作僅針對整體自由液面的波形或者固定點波高的數(shù)值結(jié)果與實驗結(jié)果的對比.而液體晃蕩最重要的危害是會產(chǎn)生瞬間增大的砰擊荷載，該荷載作用于結(jié)構(gòu)，可能會導(dǎo)致結(jié)構(gòu)失效，因此分析液體晃蕩引起的砰擊荷載更具有工程意義.Delorme等采用SPH 方法模擬了縱搖激勵作用下共振頻率附近低載液率流體晃蕩［11］，分別與實驗測量的最大沖擊壓力和沖擊壓力時程對比，但沒有考慮水平激勵下液體晃蕩荷載的特性.基于SPH 方法，Zhu等研究了在縱蕩激勵、非共振頻率和低載液率工況下，液體砰擊荷載特性［12］.

本文采用SPH 方法在二維矩形液艙內(nèi)分別針對大振幅橫蕩和橫搖激勵，研究低載液率和自由液面的一階固有頻率f1附近（共振頻率和非共振頻率）工況下液體晃蕩特性，并開展相應(yīng)的實驗研究以驗證SPH 方法模擬自由液面和砰擊荷載的合理性.

1 SPH 方法基本理論

1.1 控制方程

SPH 理論中將N－S的3個基本方程（連續(xù)方程、動量方程和狀態(tài)方程）進(jìn)行如下離散：

連續(xù)方程

式中：ρi是流體密度；Wij是與粒子i和j的位置相關(guān)的核函數(shù)；m是質(zhì)量；v是速度，vij＝vi－vj；xi表示粒子i的坐標(biāo).

動量方程

狀態(tài)方程

式（3）是動量方程的黏性形式.其中pi和pj分別是粒子i和粒子j處的壓強(qiáng)，可以通過狀態(tài)方程（5）得到；g是重力加速度.式（4）是Monaghan型人工黏性hj）；α和β是標(biāo)準(zhǔn)常數(shù)，本文中α的范圍是0.02～0.10，并且β＝0.式（5）中ρ0 為流體的初始密度；cs是聲速；γ是常數(shù)，一般取7.

1.2 邊界條件

對于固壁邊界效應(yīng)可以用多種方法來引入，本文采用動力邊界條件，式（2）、（3）和（5）這3個控制流體粒子的方程同樣可以用來描述邊界粒子，但這只是計算得到了邊界粒子合適的密度，而邊界粒子的運動僅依賴于其被定義的運動狀態(tài).為了防止邊界粒子的非物理穿透，密度低于初始密度的邊界粒子，在每一時間步中，其密度都會被自動重置成初始密度.

1.3 密度初始化

為了抑制壓力場數(shù)值計算的噪聲，引入Shepard filter算法，在每計算30 步后，通過式（6）和（7）做一次濾波處理.該算法會在當(dāng)前時間步對每個粒子做密度重分配，從而光滑壓力場.這種處理方法已經(jīng)在Colagrossi等的工作中被驗證［13］.

2 數(shù)值算例

2.1 數(shù)值模型的建立

采用SPH 方法模擬二維矩形液艙內(nèi)液體在自由液面共振頻率附近的晃蕩.數(shù)值實驗所用二維模型液艙長（L）×高（H）為97.0cm×92.7 cm.沖擊壓力監(jiān)測點P1設(shè)置在靜水面處，如圖1所示.設(shè)定流體粒子的初始間距為15 mm.計算時間步長和數(shù)據(jù)采集頻率分別為1.0×10－5s和50Hz.

2.2 實驗?zāi)Ｐ偷慕?/h3>

為了驗證SPH 方法模擬液體晃蕩的合理性，本文在大連理工大學(xué)工業(yè)裝備與結(jié)構(gòu)分析國家重點實驗室開展了一系列與數(shù)值實驗對應(yīng)的實驗研究.實驗采用的大尺度六自由度運動平臺設(shè)計承重荷載為12t，可模擬單個自由度、多個自由度耦合的規(guī)則和不規(guī)則運動.實驗系統(tǒng)如圖2所示.

圖1 模型幾何尺寸（單位：cm）Fig.1 The dimension of model（unit：cm）

圖2 晃蕩模型實驗系統(tǒng)Fig.2 The system for sloshing model experiments

為了盡可能降低數(shù)值實驗和物理實驗的差別，應(yīng)保證物理模型液艙內(nèi)流場二維性，物理模型液艙的寬度設(shè)計主要考慮以下幾個方面［6］：

（1）液艙寬度應(yīng)足夠小，以保證全局流場盡量表現(xiàn)出二維特性，避免三維流場效應(yīng)（漩渦等）的影響.

（2）液艙寬度要足夠大，以保證可以忽略沿著液艙寬度方向上產(chǎn)生的邊界層厚度對測試結(jié)果的影響［14］.

為了保證流場的二維性，歐洲多家研究機(jī)構(gòu)在全球晃蕩模型實驗基準(zhǔn)（sloshing model test benchmark，SMT）下開展了大量的數(shù)值計算，并給出SMT 實驗?zāi)Ｐ鸵号摵侠沓叽绶秶?5］.

基于上述因素，本文采用的矩形液艙寬（B）為15.8cm，物理模型的艙長和高與數(shù)值模型一致.該模型液艙采用透明的有機(jī)玻璃制成，厚度為10mm，可視為剛性模型.

2.3 模擬工況

設(shè)計的橫蕩和橫搖簡諧激勵信號如下式：

式中：ηa 為激勵振幅；f為外激激勵頻率；t為計算時間；φ為初相位；Y為激勵位移.本文中橫蕩激勵振幅η2a ＝30mm；橫搖激勵振幅η4a ＝3°，橫搖中心在艙底中心處；初相位φ＝π／2；計算時間t＝30s.

通過邊界條件、線性自由液面運動和動力條件求解二維Laplace控制方程，得到二維矩形液艙自由液面的固有頻率計算公式：

式中：g為重力加速度；h為載液高度；L為液艙的特征尺寸；n為階數(shù)；ωn為自由液面固有角頻率.

式中：fn是n階自由液面固有頻率.

靜水面高度（h）為18.6cm，即模型液艙載液率（h／H）為20%.詳細(xì)的實驗工況見表1.所有物理實驗工況均與SPH 數(shù)值實驗相同.

表1 實驗工況Tab.1 Test condition

2.4 砰擊壓力監(jiān)測

采用壓力監(jiān)測半徑內(nèi)所有粒子的壓力平均值代表監(jiān)測點的壓力值.壓力監(jiān)測半徑R為3倍的初始粒子間距（如圖3所示）.

圖3 壓力監(jiān)測點及其計算半徑Fig.3 Pressure monitoring point and its calculated radius

物理實驗中采用孔徑為5mm 的薄膜型壓阻式壓力傳感器監(jiān)測P1點處的沖擊壓力［6］.

3 結(jié)果分析與討論

3.1 自由液面波形對比

圖4給出了橫蕩激勵和1.0f1下，沖擊側(cè)壁之前的SPH 模擬的全局自由液面波形和相應(yīng)的實驗波形.在共振頻率下，液體在沖擊側(cè)壁之前，波面發(fā)生了明顯的卷波現(xiàn)象.因此，波面水頭砰擊側(cè)壁的瞬間會夾雜一些氣團(tuán).通過數(shù)值模擬發(fā)現(xiàn)氣團(tuán)僅發(fā)生在1.0f1附近.

橫搖激勵和1.1f1下，4個連續(xù)時刻中SPH模擬的全局自由液面波形和實驗結(jié)果具有較好的一致性（如圖5所示）.數(shù)值模擬的沖擊側(cè)壁前后的波面均可與實驗波形吻合.

圖5 橫搖激勵和1.1f1 下4 個連續(xù)時刻（18.0、18.2、18.4、18.6s）SPH 和 實驗的自由液面對比Fig.5 The comparison of free surface from SPH and experiments for four continuous moments （18.0，18.2，18.4，18.6 s）under 1.1f1and roll forced excitation

3.2 沖擊壓力對比

對20%的低載液率工況，液體對垂直側(cè)壁的砰擊較嚴(yán)重.分別對橫蕩和橫搖激勵下，數(shù)值計算和實驗監(jiān)測的壓力監(jiān)測點P1處的壓力時程進(jìn)行對比分析和討論.

橫蕩激勵下，4個激勵頻率在P1處監(jiān)測到的壓力時程曲線對比如圖6所示.除了共振頻率1.0f1以外，SPH 在其他激勵頻率下模擬所得壓力時程曲線與實驗結(jié)果吻合良好.在共振頻率下，數(shù)值結(jié)果與實驗結(jié)果出現(xiàn)較大差異的原因主要是因為共振頻率下，液體砰擊壁面的物理過程非常復(fù)雜，波面水頭砰擊側(cè)壁時會夾雜一些氣團(tuán)（如圖4所示）.物理實驗中，由于氣體和液體的可壓縮性，氣團(tuán)會對液體沖擊壁面起到緩沖作用，從而有效降低了沖擊壓力.而本文所采用的是單相SPH數(shù)值模擬，沒有考慮氣體相，所以液體在沖擊壁面時，就不可能有氣團(tuán)來起到緩沖作用.因此建議采用兩相流來模擬共振頻率下液體劇烈的晃蕩運動.

當(dāng)激勵頻率為1.2f1時，數(shù)值和實驗均體現(xiàn)出“拍”現(xiàn)象（beating effect）.此外，在非共振頻率1.1f1處，放大一個沖擊過程，如圖7所示.SPH可以模擬砰擊荷載的雙峰沖擊特性.第一個峰值代表液體沖擊艙壁產(chǎn)生的砰擊荷載（動壓力）.第二個峰值是由于液體砰擊壁面后，沿著艙壁爬升，當(dāng)運動速度為零時，動能最小，勢能最大導(dǎo)致的.因此第二個峰值表征靜水力.此外，從圖7中可以看出，在同一個周期內(nèi)，SPH 數(shù)值模擬得到的壓力第二個峰值要高于實驗結(jié)果，而且這種現(xiàn)象普遍存在于相應(yīng)激勵頻率下的整個沖擊過程.出現(xiàn)這種現(xiàn)象的主要原因，是由于實驗所用的模型液艙具有一定寬度，因此會有三維效應(yīng).通過實驗觀測發(fā)現(xiàn)液體在沿垂直壁面爬升運動過程中會出現(xiàn)沿寬度方向的發(fā)散運動，這種發(fā)散會使得回落到P1位置處的液體水團(tuán)質(zhì)量減小，從而減低了液體水團(tuán)回落時對測點P1處的沖擊壓力.

圖7 橫蕩激勵和1.1f1 下壓力監(jiān)測點P1處的壓力時程放大圖Fig.7 Enlarged view of the pressure time history at P1under 1.1f1and sway forced excitation

橫搖激勵下，壓力監(jiān)測點P1 處的壓力時程對比如圖8所示.橫搖激勵的沖擊荷載較橫蕩激勵下更顯著.SPH 模擬的數(shù)值結(jié)果較實驗結(jié)果要大，分析認(rèn)為，主要原因是液艙三維效應(yīng)及液體劇烈砰擊艙壁時瞬間產(chǎn)生的氣體和液體壓縮作用.

4 結(jié) 論

（1）SPH 方法可以很好地模擬低載液率橫蕩和橫搖激勵下全局自由液面的波形，如水躍、破碎波等自由液面的大變形現(xiàn)象.

（2）SPH 方法模擬共振頻率下的砰擊荷載比實驗結(jié)果要大，主要原因是本文忽略氣體相影響.因此建議采用兩相流來模擬共振頻率下液體劇烈的晃蕩運動.

（3）SPH 方法可以很好地模擬橫蕩和橫搖激勵下非共振頻率處晃蕩荷載的主要特點.因此，該方法可用于模擬和評估非共振區(qū)域晃蕩荷載的特性.

（4）采用SPH 方法對共振頻率激勵下儲液艙內(nèi)液體晃蕩荷載的特性進(jìn)行研究發(fā)現(xiàn)，液體砰擊艙壁前會發(fā)生水躍現(xiàn)象，且晃蕩沖擊荷載具有雙峰沖擊特性.

［1］ 衛(wèi)志軍，岳前進(jìn)，阮詩倫，等.矩形液艙晃蕩沖擊載荷的試驗機(jī)理研究［J］.船舶力學(xué)，2012，16（8）：885－892.WEI Zhi－jun，YUE Qian－jin，RUAN Shi－lun，etal.An experimental investigation of liquid sloshing impact load on a rectangular tank［J］.Journal of Ship Mechanics，2012，16 （8）：885－892.（in Chinese）

［2］ ZHU Ren－qing，WU You－sheng，Incecik A.Numerical simulation of liquid sloshing—a review［J］.Shipbuilding of China，2004，45（2）：14－27.

［3］ Lucy L B.A numerical approach to the testing of the fission hypothesis ［J］.The Astronomical Journal，1977，82：1013－1024.

［4］ Gingold R A，Monaghan J J.Smoothed particle hydrodynamics— Theory and application to nonspherical stars［J］.Monthly Notices of the Royal Astronomical Society，1977，181：375－389.

［5］ Monaghan J J.Simulating free surface flows with SPH ［J］.Journal of Computational Physics，1994，110（2）：399－406.

［6］ 衛(wèi)志軍，岳前進(jìn)，張文首，等.大尺度儲艙液體晃蕩砰擊壓力測量方法研究［J］.中國科學(xué)：物理學(xué) 力學(xué) 天文學(xué)，2014，44（7）：746－758.WEI Zhi－jun，YUE Qian－jin，ZHANG Wen－shou，etal.Experimental investigation of violent slamming pressure in large－scaled tank［J］.Scientia Sinica：Physica，Mechanica ＆ Astronomica，2014，44（7）：746－758.（in Chinese）

［7］ 陳正云，朱仁慶，祁江濤.基于SPH 法的二維液體大幅晃蕩數(shù)值模擬［J］.船海工程，2008，37（2）：44－47.CHEN Zheng－yun，ZHU Ren－qing，QI Jiang－tao.Numerical simulation of sloshing in two dimensional liquid tank based on SPH method ［J］.Ship ＆Ocean Engineering，2008，37（2）：44－47.（in Chinese）

［8］ 崔 巖，吳 衛(wèi)，龔 凱，等.二維矩形水槽晃蕩過程的SPH 方法模擬［J］.水動力學(xué)研究與進(jìn)展，2008，23（6）：618－624.CUI Yan，WU Wei，GONG Kai，etal.Numerical simulation of sloshing in two dimensional rectangular tanks with SPH ［J］.Chinese Journal of Hydrodynamics，2008，23 （6）：618－624.（in Chinese）

［9］ 李大鳴，陳海舟，張建偉，等.基于SPH 法的二維矩形液艙晃蕩研究［J］.計算力學(xué)學(xué)報，2010，27（2）：369－374.LI Da－ming，CHEN Hai－zhou，ZHANG Jian－wei，etal.A study of 2Dliquid sloshing in rectangle tanks based on SPH method［J］.Chinese Journal of Computational Mechanics，2010，27（2）：369－374.（in Chinese）

［10］ 李大鳴，李玲玲，高 偉，等.液艙內(nèi)液體晃蕩的SPH 模型［J］.哈爾濱工程大學(xué)學(xué)報，2012，33（1）：37－41.LI Da－ming，LI Ling－ling，GAO Wei，etal.A model of liquid sloshing in tanks based on the smoothed particle hydrodynamics （SPH）method［J］.Journal of Harbin Engineering University，2012，33（1）：37－41.（in Chinese）

［11］ Delorme L，Colagrossi A，Souto－Iglesias A，etal.A set of canonical problems in sloshing，Part I：Pressure field in forced roll—comparison between experimental results and SPH ［J］.Ocean Engineering，2009，36（2）：168－178.

［12］ ZHU Ren－qing，CHEN Zheng－yun，WANG Quan.Numerical simulation of 2Dsloshing in liquid tanks based on SPH method ［C］／／ASME 2010 29th International Conference on Ocean，Offshore and Arctic Engineering，OMAE2010.New York：American Society of Mechanical Engineers，2010：961－966.

［13］ Colagrossi A，Landrini M.Numerical simulation of interfacial flows by smoothed particle hydrodynamics ［J］.Journal of Computational Physics，2003，191（2）：448－475.

［14］ Faltinsen O M，F(xiàn)iroozkoohi R，Timokha A N.Steady－state liquid sloshing in a rectangular tank with a slat－type screen in the middle：Quasilinear modal analysis and experiments ［J］.Physics of Fluids，2011，23（4）：042101.

［15］ Loysel T，Chollet S，Gervaise E，etal.Results of the first sloshing model test benchmark ［C］／／Proceedings of the 22nd （2012）International Offshore and Polar Engineering Conference.Mountain View：International Society of Offshore and Polar Engineers，2012：398－408.