SPH求解器的開發(fā)及其在液艙晃蕩數(shù)值模擬中的應(yīng)用

尹梓煒 陳曉明 趙成璧

摘 要：本文提出一種改良的SPH固壁邊界條件處理方法，并基于此開發(fā)適用于模擬液艙晃蕩的數(shù)值求解器。采取該求解器對(duì)液艙晃蕩現(xiàn)象進(jìn)行數(shù)值模擬，結(jié)果與相關(guān)文獻(xiàn)結(jié)果相當(dāng)吻合，驗(yàn)證了其可靠性和準(zhǔn)確性。在此基礎(chǔ)上，通過計(jì)算分析橫蕩狀態(tài)下不同隔板高度液艙晃蕩時(shí)所產(chǎn)生的壁面壓力和液體重心高度變化趨勢(shì)，探討了隔板高度對(duì)液艙晃蕩的影響規(guī)律和物理機(jī)制，并為防蕩隔板的設(shè)計(jì)提出合理建議。

關(guān)鍵詞：SPH；固壁邊界條件；液艙晃蕩；數(shù)值模擬

中圖分類號(hào)：U661.74 文獻(xiàn)標(biāo)識(shí)碼：A

Abstract： A modified treatment for rigid boundary condition of SPH method is proposed in this paper， and a numerical solver that suitable for sloshing simulation is developed based on this SPH method. The phenomenon of liquid sloshing in a tank is simulated by using this solver， and the numerical simulation results are in good agreement with the results from the related literature， which verifies the reliability and accuracy of the model. On the basis of that， through the calculation and analysis of the wall pressure and the trend of liquid gravity center height， which generated by sloshing with different baffle lengths under the sway condition， this paper preliminarily discusses the effect and physical mechanism of baffle length on sloshing and proposes reasonable suggestions about the design of anti-sloshing baffle.

Key words： SPH； Rigid boundary condition； Sloshing； Numerical simulation

1 引言

液艙晃蕩是指在部分充液的容器中，其內(nèi)部液體在外部激勵(lì)下產(chǎn)生運(yùn)動(dòng)的現(xiàn)象。液艙晃蕩是船舶航行中遇到的常見現(xiàn)象，具有較強(qiáng)的隨機(jī)性和非線性。液體晃蕩產(chǎn)生的瞬時(shí)壓力可對(duì)艙壁結(jié)構(gòu)造成損傷，同時(shí)晃蕩產(chǎn)生的附加力矩會(huì)影響船舶穩(wěn)性，危害船舶安全，因此液艙防晃是船舶工程領(lǐng)域的基礎(chǔ)性課題，具有重要的工程應(yīng)用價(jià)值。

液艙晃蕩問題的研究方法主要分為實(shí)驗(yàn)、理論分析和數(shù)值仿真。隨著計(jì)算機(jī)科學(xué)的發(fā)展，數(shù)值分析已經(jīng)成為研究液艙晃蕩問題的重要手段。傳統(tǒng)的數(shù)值方法如MAC、VOF等是基于網(wǎng)格對(duì)流場(chǎng)進(jìn)行求解，在處理液艙晃蕩這類帶自由液面的強(qiáng)非線性運(yùn)動(dòng)問題時(shí)存在困難。光滑粒子流體動(dòng)力學(xué)方法（SPH），最早由Lucy[1]和Monaghan[2]等分別獨(dú)立提出用于解決天體物理學(xué)問題，后經(jīng)過改進(jìn)引入到水動(dòng)力學(xué)中。SPH作為一種無(wú)網(wǎng)格的拉格朗日方法，它克服了無(wú)網(wǎng)格方法在處理流體自由面大變形問題上的弱點(diǎn)，近年來(lái)逐漸被應(yīng)用于液艙晃蕩問題的研究。

較早期的有Iglesias[3][4]等利用SPH方法先后研究了減搖水艙晃蕩特性和不同頻率下液艙晃蕩力矩幅值等問題，后來(lái)崔巖[5]等運(yùn)用SPH方法模擬分析了二維矩形水槽在激勵(lì)頻率接近一階固有頻率的縱蕩過程中的晃蕩現(xiàn)象。在基于SPH方法的液艙晃蕩數(shù)值模擬中，邊界條件的處理十分關(guān)鍵，而上述文獻(xiàn)均采用Lennard-Jones斥力模型實(shí)施固壁邊界條件。該模型對(duì)靠近邊界的粒子施加強(qiáng)斥力，這容易導(dǎo)致粒子的初始分布被破壞，且時(shí)間步長(zhǎng)受到嚴(yán)格的限制。同時(shí)為了防止粒子穿透邊界，特別是在邊界形狀突變處，往往使用多層邊界粒子，降低了計(jì)算效率，且會(huì)導(dǎo)致邊界附近產(chǎn)生明顯的壓力振蕩。為改善該問題，Shao[6]等提出耦合動(dòng)力學(xué)邊界的SPH方法并研究了二維液艙晃蕩問題，得到較為精確的流場(chǎng)壓力。Chen[7]等在此基礎(chǔ)上修正其人工斥力的表達(dá)式，模擬了二維矩形容器在液艙晃蕩下的艙壁砰擊壓力并進(jìn)行了實(shí)驗(yàn)對(duì)比，獲得了讓人滿意的結(jié)果。

基于上述研究成果，為了獲得更為精確的液艙晃蕩模擬結(jié)果，本文在提出一種改良的固壁邊界條件處理方法基礎(chǔ)上，自主開發(fā)了適用于液艙晃蕩的SPH數(shù)值求解器，并進(jìn)一步研究了隔板高度對(duì)晃蕩行為的影響，從而擴(kuò)展了SPH方法在液艙晃蕩中的應(yīng)用。

2 SPH求解器的開發(fā)

2.1 SPH基本方法和理論

SPH方法中[8]，計(jì)算域被離散為一系列具有相互作用的粒子，某時(shí)刻粒子的物理量可由一定光滑半徑內(nèi)的周圍粒子通過核函數(shù)（如三次樣條函數(shù)等）估計(jì)得到?；赟PH的核近似和粒子近似思想，考慮重力下控制流體運(yùn)動(dòng)的Navier-Stokes方程可離散為如下形式：

壓力通過下列狀態(tài)方程求得：

2.2 自適應(yīng)邊界粒子法向斥力模型

Akinic[10]提出了一種流固耦合邊界方法，該方法引入邊界粒子質(zhì)量函數(shù)，考慮邊界粒子的影響對(duì)流體粒子的密度和受力進(jìn)行修正。本文基于該方法的思想并對(duì)其進(jìn)行改進(jìn)，只考慮法向斥力，提出一種改良的自適應(yīng)邊界粒子法向斥力模型，得到考慮邊界粒子b的貢獻(xiàn)的SPH控制方程如下：

由上式可知，由于核函數(shù)隨距離增加而單調(diào)遞減的性質(zhì)，邊界粒子分布密集的地方不會(huì)產(chǎn)生過大斥力，分布稀疏的地方也不會(huì)產(chǎn)生過小的斥力，因此相對(duì)于傳統(tǒng)SPH邊界處理方法具有自適應(yīng)性。且該方法斥力大小又與流體粒子自身壓力大小有關(guān)，只在法向施加排斥力，減少邊界粒子對(duì)流場(chǎng)的擾動(dòng)，可以使粒子在邊界處受力更均勻。

以簡(jiǎn)單的潰壩模型為算例。從圖1的速度矢量對(duì)比可以看出：傳統(tǒng)的Lennard-Jones斥力模型中潰壩水頭前沿的粒子運(yùn)動(dòng)十分不均勻，明顯是由于產(chǎn)生的瞬時(shí)斥力過大導(dǎo)致；而本文方法中的底部邊界粒子產(chǎn)生足夠而均勻的斥力，有效避免粒子在固壁邊界發(fā)生非物理穿透的現(xiàn)象，同時(shí)使粒子運(yùn)動(dòng)平滑有序，證明本文方法的效果明顯優(yōu)于傳統(tǒng)的斥力模型。

2.3 算法流程

SPH求解器的算法流程如圖2。

其中：鄰域粒子搜索方法有全配對(duì)搜索法、鏈表搜索法和樹形搜索法等方法；時(shí)間積分的方法可采用蛙跳法、預(yù)測(cè)校正法和龍格庫(kù)塔法等方法。本文采用鏈表搜索法進(jìn)行鄰域搜索，采用蛙跳法進(jìn)行時(shí)間積分，時(shí)間步長(zhǎng)滿足CFL條件，具體見參考文獻(xiàn)[8]。

3 基于SPH求解器的液艙晃蕩數(shù)值模擬

3.1 典型液艙晃蕩實(shí)驗(yàn)

Chen[7]的實(shí)驗(yàn)是研究液艙晃蕩的典型實(shí)驗(yàn)，并給出了較準(zhǔn)確的SPH模擬結(jié)果。為了驗(yàn)證本文液艙晃蕩SPH方法數(shù)值求解器的可靠性和精確性，本文模擬了Chen液艙晃蕩算例，并與其實(shí)驗(yàn)結(jié)果和模擬結(jié)果進(jìn)行比較。

計(jì)算模型中矩形液艙長(zhǎng)和高為L(zhǎng)=H=1 m、充液深度D=0.3 m、橫搖頻率3.81rad/s、振幅5°，容器底部中心處為軸心。計(jì)算采用流體粒子數(shù)12 870、時(shí)間步長(zhǎng)1×10-4 s進(jìn)行模擬，流體壓縮率控制在0.1%附近。

圖3和圖4分別為矩形液艙在激勵(lì)頻率3.81 rad/s的橫搖下晃蕩的四個(gè)典型時(shí)刻的chen實(shí)驗(yàn)結(jié)果和本文模擬結(jié)果。從圖中可以看出，本文的SPH模型有效地模擬了自由液面形態(tài)，同時(shí)獲得了較精確的均勻壓力分布。

在本文中，基于SPH的核近似和粒子近似的思想，邊界粒子b的壓強(qiáng)通過周圍流體粒子根據(jù)

插值估算。取距底部0.2 m高度處艙壁作為壓力監(jiān)測(cè)點(diǎn)。

圖5所示為壓力監(jiān)測(cè)點(diǎn)在計(jì)算時(shí)段2～10 s的壓力時(shí)間歷程。由圖5可知，本文SPH方法的壓力計(jì)算結(jié)果與Chen的數(shù)值模擬結(jié)果的變化趨勢(shì)基本一致，同時(shí)壓力的峰值則與Chen的實(shí)驗(yàn)結(jié)果更為接近，證明本文SPH方法在液艙晃蕩的數(shù)值模擬應(yīng)用具有更高的可靠性和精確性。

3.2 防蕩隔板高度對(duì)液艙晃蕩影響

為進(jìn)一步研究不同液艙結(jié)構(gòu)對(duì)晃蕩行為的影響，本文在與3.1算例尺寸相同的矩形液艙中添加防蕩隔板，模擬了不同高度的I型單隔板在橫蕩頻率接近液艙一階模態(tài)的橫蕩工況下的防蕩效果對(duì)比。

根據(jù)線性理論，由 可得矩形液艙晃蕩的一階固有頻率為ω1=4.76 rad/s，對(duì)液艙施加簡(jiǎn)諧激勵(lì)x=Asin（ωt），取ω=1.166ω1、振幅A=0.05 m。在橫蕩工況下，為達(dá)到防蕩效果最大化，在液艙中央豎向布置I型隔板，隔板高度l分別取0.25 D、0.5 D和0.75 D。同時(shí)為了分析方便，取距離艙壁底部0.3 m處側(cè)壁作為壓力監(jiān)測(cè)點(diǎn)。

圖6（1）描述了不同單隔板尺寸工況下艙壁壓力計(jì)算監(jiān)測(cè)點(diǎn)在計(jì)算時(shí)段2～10 s的壓力時(shí)間歷程。圖6（2）描述了不同單隔板尺寸工況下在計(jì)算時(shí)段0～10 s液體重心高度的時(shí)間歷程。由圖6（1）可知，隨著隔板高度的增加，整體上壓力峰值呈減少趨勢(shì)且雙峰特性減弱，當(dāng)隔板高度為0.75 D時(shí)最大壓力峰值能降至無(wú)隔板狀態(tài)下的24.5%。同時(shí)，當(dāng)隔板增加到一定高度后，壓力峰值對(duì)隔板高度變化的敏感度降低。由圖6（2）可知，隔板高度和液體重心高度呈明顯的正相關(guān)，其中當(dāng)隔板高度為0.75D時(shí)，重心高度波動(dòng)最大幅值可降至無(wú)隔板狀態(tài)時(shí)的19.2%。

為了進(jìn)一步研究液艙晃蕩時(shí)隔板高度對(duì)壓力峰值影響的物理機(jī)理，我們對(duì)不同隔板高度下的液艙流場(chǎng)速度進(jìn)行對(duì)比分析。

圖7反映了t=2.26 s時(shí)刻不同隔板高度的流場(chǎng)速度矢量和艙壁壓力?？梢钥闯觯弘S著隔板高度的增加，隔板附近流場(chǎng)受到干擾，出現(xiàn)了紊亂并逐漸產(chǎn)生了渦。此時(shí)，渦街阻尼效應(yīng)對(duì)較大波高的形成產(chǎn)生了削弱作用，有效抑制了液艙晃蕩時(shí)的液面劇烈波動(dòng)，從而有效地降低艙壁晃蕩時(shí)的艙壁壓力。

不過，結(jié)合圖6（1）可知，當(dāng)隔板到達(dá)一定高度后，其對(duì)艙壁壓力峰值的影響效果將減弱；同時(shí)圖7表明，隨著隔板高度的增加，隔板的渦街振動(dòng)響應(yīng)不斷增強(qiáng)，隔板的液體沖擊荷載作用不斷增大，這會(huì)導(dǎo)致隔板端部應(yīng)力集中，因此從安全性和經(jīng)濟(jì)型考慮，隔板高度應(yīng)適中，建議隔板在具體設(shè)計(jì)時(shí)可采用本文數(shù)值模擬方法確定隔板高度尺寸。

4 結(jié)論

本文利用自主開發(fā)改進(jìn)的SPH求解器研究了二維矩形液艙晃蕩問題，取得較好的應(yīng)用效果。由數(shù)值模擬結(jié)果可以看出：

（1）本文改良的固壁邊界條件處理方法，能有效阻止粒子穿透邊界并使粒子在邊界處的運(yùn)動(dòng)均勻有序。以此為基礎(chǔ)自主開發(fā)的SPH求解器具有較高的可靠性和精確性，能很好地捕捉液艙晃蕩時(shí)的自由液面狀態(tài)，獲得較為準(zhǔn)確的流場(chǎng)壓力值，是對(duì)液艙晃蕩現(xiàn)象有效的數(shù)值模擬手段，具有較好的應(yīng)用前景；

（2）本文進(jìn)一步的研究結(jié)果表明，防蕩隔板能明顯降低液艙晃蕩的幅度和艙壁壓力峰值，但從安全性和經(jīng)濟(jì)性的角度出發(fā)，隔板高度應(yīng)適中，不宜過高。

此外，SPH方法中在增加算法穩(wěn)定性、提高計(jì)算精度和效率等方面尚有許多需要完善之處，本文的研究為進(jìn)一步開發(fā)高效通用的SPH求解器以及拓展其應(yīng)用領(lǐng)域打下了基礎(chǔ)。

參考文獻(xiàn)

[1]Lucy L B. A numerical approach to the testing of the fission hypothesis[J]. Astronomical Journal， 1977， （82）.

[2]Gingold R A， Monaghan J J. Smoothed particle hydrodynamics： theory and application to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society， 1977， 181（3）.

[3]Souto Iglesias A， Pérez Rojas L， Zamora Rodríguez R. Simulation of anti-roll tanks and sloshing type problems with smoothed particle hydrodynamics[J]. Ocean Engineering， 2004， 31（8）.

[4]Souto-Iglesias A， Delorme L， Pérez-Rojas L， et al. Liquid moment amplitude assessment in sloshing type problems with smooth particle hydrodynamics[J]. Ocean Engineering， 2006， 33（11）.

[5]崔巖， 吳衛(wèi)， 劉樺. SPH方法模擬二維矩形水槽晃蕩過程[C]// 全國(guó)水動(dòng)力學(xué)研討會(huì). 2006.

[6]Shao J R， Li H Q， Liu G R， et al. An improved SPH method for modeling liquid sloshing dynamics[J]. Computers & Structures， 2012， s 100–101（6）.

[7]Chen Z， Zong Z， Li H T， et al. An investigation into the pressure on solid walls in 2D sloshing using SPH method[J]. Ocean Engineering， 2013， 59（2）.

[8]G.R.Liu， M.B.Liu， Liu，等. 光滑粒子流體動(dòng)力學(xué)：一種無(wú)網(wǎng)格粒子法[M]. 湖南大學(xué)出版社， 2005.

[9]Monaghan J J， Gingold R A. Shock simulation by the particle method SPH[J]. Journal of Computational Physics， 1983， 52（2）.

[10]Akinci N， Ihmsen M， Akinci G， et al. Versatile rigid-fluid coupling for incompressible SPH[J]. Acm Transactions on Graphics， 2012， 31（4）.