動(dòng)態(tài)邊界粒子SPH法在養(yǎng)殖工船橫搖液艙晃蕩模擬中的應(yīng)用比較

趙新穎，黃溫赟，黃文超，管延敏

(1 中國(guó)水產(chǎn)科學(xué)研究院漁業(yè)機(jī)械儀器研究所，農(nóng)業(yè)農(nóng)村部遠(yuǎn)洋漁船與裝備重點(diǎn)實(shí)驗(yàn)室，上海 200092;2 青島海洋科學(xué)與技術(shù)國(guó)家試點(diǎn)實(shí)驗(yàn)室深藍(lán)漁業(yè)工程聯(lián)合實(shí)驗(yàn)室，山東 青島 266237;3 江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003)

近年來(lái)，作為漁業(yè)轉(zhuǎn)型升級(jí)和綠色發(fā)展的重要方向，工船養(yǎng)殖受到廣泛關(guān)注，養(yǎng)殖水艙是養(yǎng)殖工船的核心部分，與載液船舶相似，養(yǎng)殖水艙在波浪作用下產(chǎn)生的液艙晃蕩現(xiàn)象，對(duì)船體性能以及液艙艙壁結(jié)構(gòu)強(qiáng)度產(chǎn)生影響，進(jìn)而影響?zhàn)B殖安全。為了對(duì)此類現(xiàn)象進(jìn)行研究，有限差分法[1]、有限體積法[2]、有限元法[3]、邊界元法[4]、無(wú)網(wǎng)格法(MPS[5]、SPH[6]、CPM[7])等諸多方法相繼被用于晃蕩現(xiàn)象數(shù)值預(yù)報(bào)。劇烈的液體晃蕩會(huì)產(chǎn)生的自由液面翻卷、破碎等現(xiàn)象的模擬給傳統(tǒng)有網(wǎng)格方法帶來(lái)一定的難度，光滑粒子流體動(dòng)力學(xué)(SPH)法作為無(wú)網(wǎng)格法的一種，擺脫了傳統(tǒng)方法對(duì)網(wǎng)格、單元的依賴，避免了拉格朗日方法中的網(wǎng)格扭曲及網(wǎng)格重構(gòu)問(wèn)題。該方法通過(guò)粒子的運(yùn)動(dòng)模擬流場(chǎng)流動(dòng)現(xiàn)象，可以實(shí)現(xiàn)對(duì)自由液面的精確追蹤[1]，在處理大變形[9]、動(dòng)邊界[10]、自由液面翻卷、破碎[11-12]等問(wèn)題時(shí)具有明顯的優(yōu)勢(shì)，因此SPH已成為研究液艙晃蕩問(wèn)題的熱門方法。如李大鳴等[13]對(duì)中、高液位液體晃蕩及防蕩隔板對(duì)液體晃蕩現(xiàn)象的抑制作用進(jìn)行了數(shù)值模擬，Shao等[14]提出了一種改進(jìn)SPH方法用于模擬液體晃蕩現(xiàn)象，Gotoh等[15]運(yùn)用不可壓縮SPH法模擬了劇烈液體晃蕩現(xiàn)象，曹雪雁等[16-17]運(yùn)用SPH法模擬了矩形箱體液體晃蕩現(xiàn)象以及船舶破艙進(jìn)水特性。

SPH方法中，壁面邊界的處理是提高SPH的精度和穩(wěn)定性的重要因素，一直是SPH研究關(guān)注的熱點(diǎn)[18]。位于固壁邊界上或邊界附近的粒子由于其核函數(shù)被邊界截?cái)?，使得該區(qū)域的變量梯度計(jì)算不準(zhǔn)確，產(chǎn)生截?cái)嗾`差從而引起數(shù)值結(jié)果畸變，無(wú)法保證二階精度[19]。近年來(lái)，不同的壁面邊界的處理方法不斷被提出，傳統(tǒng)的鏡像對(duì)稱邊界法[20]受鏡像對(duì)稱算法的限制，一般只應(yīng)用于規(guī)則的平面或者直角邊界，在復(fù)雜邊界問(wèn)題中很少被使用；Monaghan[21]提出的排斥力法方法簡(jiǎn)單，但因守恒性差，無(wú)法精確模擬粒子在邊界附近的運(yùn)動(dòng)，影響了壁面載荷計(jì)算精度；Libersky等[22]引入的鏡像粒子法能夠?qū)?nèi)部粒子進(jìn)行準(zhǔn)確計(jì)算，但可能會(huì)產(chǎn)生局部流體粒子穿透邊界現(xiàn)象；Marrone等[12]在鏡像粒子的基礎(chǔ)上提出的固定虛粒子法，滿足了壁面不可穿透條件，在處理復(fù)雜邊界問(wèn)題時(shí)，依然存在虛粒子布置的困難；動(dòng)態(tài)邊界粒子由Crespo等[24]詳細(xì)介紹，并在Liu等[25]的研究中被進(jìn)一步使用，該方法虛粒子滿足與流體粒子相同的連續(xù)性和狀態(tài)方程，可以與流體粒子在同一循環(huán)內(nèi)簡(jiǎn)單地進(jìn)行計(jì)算，可以節(jié)省大量的計(jì)算時(shí)間。

為準(zhǔn)確模擬養(yǎng)殖工船養(yǎng)殖艙自由液面變化以及計(jì)算艙壁受力情況，保證養(yǎng)殖工船結(jié)構(gòu)強(qiáng)度和養(yǎng)殖安全，本研究采用動(dòng)態(tài)邊界粒子SPH法對(duì)養(yǎng)殖艙橫搖激勵(lì)下液體晃蕩現(xiàn)象展開(kāi)模擬，并對(duì)動(dòng)態(tài)邊界粒子處理的兩種方法展開(kāi)比較研究，數(shù)值結(jié)果表明本文采取的方法精確有效，可以為養(yǎng)殖水艙晃蕩現(xiàn)象模擬提供數(shù)值參考。

1 SPH法數(shù)值模型

1.1 控制方程及其離散形式

對(duì)于拉格朗日形式下的弱可壓縮流體，其連續(xù)性方程和動(dòng)量方程可以表示為：

(1)

(2)

式中：p、ρ、u、μ分別表示流體的壓力、密度、速度和黏性；g為重力加速度。流體的弱可壓縮性采用Monaghan等[21]提出的人工壓縮法，如式(3)。

(3)

利用SPH法的基本原理，連續(xù)性方程和動(dòng)量方程可以離散為：

(4)

(5)

式中：pi、ρi、ui分別是第i個(gè)粒子的壓力、密度和速度；mj、ρj為粒子i支持域中相鄰粒子j的質(zhì)量和密度。ij為人工黏性項(xiàng)[26]，其表達(dá)式為：

(6)

1.2 核函數(shù)

核函數(shù)對(duì)數(shù)值模擬的精度、效率及穩(wěn)定性，目前應(yīng)用較為廣泛的核函數(shù)有B-樣條核函數(shù)、3次樣條核函數(shù)、5次樣條核函數(shù)和高斯核函數(shù)等[27]。Dehnen等[28]研究發(fā)現(xiàn)，隨著光滑長(zhǎng)度和周圍粒子數(shù)量的增加，傳統(tǒng)的核函數(shù)受到聚集不穩(wěn)定因素粒子容易產(chǎn)生結(jié)對(duì)現(xiàn)象，推薦了Wendland核函數(shù)[29]。本研究在計(jì)算過(guò)程中采用Wendland的C2核函數(shù)：

(7)

式中；對(duì)于二維問(wèn)題αd=7/4πh2。

1.3 時(shí)間積分

本研究采用預(yù)估校正法求解控制方程，計(jì)算粒子位置、密度和速度。

在預(yù)估步有：

(8)

在校正步有：

(9)

(10)

2 壁面邊界處理的動(dòng)態(tài)邊界條件

2.1 動(dòng)態(tài)邊界條件基本原理

動(dòng)態(tài)邊界條件是在邊界處布置一組動(dòng)態(tài)邊界粒子(圖1)，這些邊界粒子跟流體粒子一樣參與連續(xù)性方程計(jì)算，不參與動(dòng)量方程求解，當(dāng)流體粒子靠近邊界與邊界粒子的距離減小至2倍光滑長(zhǎng)度以下時(shí)，受流體粒子影響，邊界粒子密度增加導(dǎo)致壓力增大，邊界粒子壓力的變化又反作用于流體粒子動(dòng)量方程，引起流體粒子速度的變化，從而防止流體粒子穿透邊界。動(dòng)態(tài)邊界粒子位置可以保持固定不動(dòng)，也可以通過(guò)某種方式使其運(yùn)動(dòng)，本研究針對(duì)這兩種方式展開(kāi)比較分析。

圖1 動(dòng)態(tài)邊界條件示意圖Fig.1 Schematic diagram of dynamic particle condition

2.2 動(dòng)態(tài)邊界粒子運(yùn)動(dòng)

假定壁面平移速度為u0，繞軸心O旋轉(zhuǎn)的角速度為ω，對(duì)于任一邊界粒子，其運(yùn)動(dòng)速度uw為

uw=u0+ω×(rw-rO)

(11)

式中：rw、rO分別為邊界粒子、旋轉(zhuǎn)軸心的空間位置。

在預(yù)估步，邊界粒子的位置為：

(12)

在校正步，邊界粒子的位置為：

(13)

在每一預(yù)估步及校正步，先進(jìn)行邊界粒子速度、位置的更新，再求解流體粒子連續(xù)性方程、動(dòng)量方程。

2.3 動(dòng)態(tài)邊界粒子固定不動(dòng)

以物體壁面隨體坐標(biāo)系為參考坐標(biāo)系，保持壁面邊界不動(dòng)，流場(chǎng)所受重力方向發(fā)生變化，對(duì)于二維問(wèn)題，有如下公式：

(14)

式中：gx、gy為重力g的兩個(gè)分量；θ為物體繞軸心O旋轉(zhuǎn)的角位移。

任一流場(chǎng)粒子a所受質(zhì)量力為：

(15)

式中：第一項(xiàng)為角加速度項(xiàng)，第二項(xiàng)為向心加速度項(xiàng)，第三項(xiàng)為科氏加速度項(xiàng)。

將公式(15)代入公式(5)，動(dòng)量方程變?yōu)椋?/p>

(16)

3 數(shù)值算例

3.1 輕微液體晃蕩算例分析

通常來(lái)講，水平激勵(lì)的幅值小于15%的箱體內(nèi)液面半徑時(shí)為小振幅激勵(lì)，即輕微晃蕩；反之，則為大振幅激勵(lì)，相應(yīng)地為劇烈晃蕩。針對(duì)輕微液體晃蕩，分別采用兩種動(dòng)態(tài)邊界粒子模型對(duì)其進(jìn)行數(shù)值模擬。計(jì)算模型定義如圖2所示，矩形箱體的長(zhǎng)度L=0.92 m，高度H=0.62 m，水深h=0.465 m。不考慮空氣影響，水箱里的水密度為ρ=1 000 kg/m3。對(duì)箱體施加橫搖激勵(lì)使其繞箱體中心(軸x=L/2,y=H/2)轉(zhuǎn)動(dòng)，橫搖角位移θr=θ0sinωrt，式中，θ0為最大角位移，ωr為橫搖圓頻率。

圖2 矩形箱體坐標(biāo)系設(shè)置示意圖Fig.2 Schematic diagram of sloshing tank and the co-ordinate system

在水箱右壁面距離底部0.17 m處設(shè)置一個(gè)壓力監(jiān)測(cè)點(diǎn)A(x=0.92 m,y=0.17 m)，并通過(guò)與文獻(xiàn)[31]的實(shí)驗(yàn)值進(jìn)行比較對(duì)該算例進(jìn)行驗(yàn)證。圖3為在θ0=4°、ωr=2.0 rad/s工況下，分別采用兩種動(dòng)態(tài)邊界粒子運(yùn)動(dòng)方法A點(diǎn)壓力隨時(shí)間的變化曲線以及與實(shí)驗(yàn)結(jié)果的對(duì)比。

圖3 點(diǎn)A處壓強(qiáng)比較Fig.3 Pressure comparison at point A

由圖3可見(jiàn)，兩種方法得出的計(jì)算結(jié)果趨勢(shì)基本一致，但相對(duì)于邊界粒子運(yùn)動(dòng)法，邊界粒子固定不動(dòng)法計(jì)算值與文獻(xiàn)[31]更為接近。

為了進(jìn)一步比較兩種方法對(duì)計(jì)算精度的影響，對(duì)作用在箱體上的水動(dòng)力和旋轉(zhuǎn)繞軸的力矩進(jìn)行了對(duì)比，如圖4所示。由圖4可見(jiàn)，邊界粒子運(yùn)動(dòng)法在初始計(jì)算時(shí)會(huì)產(chǎn)生較大數(shù)值振蕩，在t=1.0 s左右計(jì)算趨于穩(wěn)定，兩種計(jì)算方法得到的水動(dòng)力和力矩的變化趨勢(shì)基本一致，兩種方法計(jì)算值大小略有差異。

圖4 作用于箱體上的力和力矩Fig.4 Time history of hydrodynamic force and moment

圖5為兩種方法在模擬過(guò)程中幾個(gè)典型時(shí)刻壓力場(chǎng)的對(duì)比，左邊為邊界粒子運(yùn)動(dòng)的壓力場(chǎng)，右邊為邊界粒子不動(dòng)的壓力場(chǎng)?？傮w來(lái)看兩種方法的模擬結(jié)果接近，均能很好地捕捉自由液面形狀。

圖5 計(jì)算結(jié)果比較Fig.5 Comparison of the numerical results between moving boundary particles(left)and fixed particles(right)

3.2 劇烈液體晃蕩算例分析

定義如圖6所示的矩形箱體，箱體寬L=0.90 m，高H=0.508 m，h=0.093 m對(duì)箱體施加橫搖激勵(lì)使其繞箱體中心(軸x=0,y=0)做簡(jiǎn)諧振動(dòng)，根據(jù)箱體自振頻率表達(dá)式：

圖6 矩形箱體坐標(biāo)系設(shè)置示意圖Fig.6 Schematic diagram of sloshing tank and the co-ordinate system

(17)

得ω0=1.919 rad/s，選取ωr=1.631 rad/s(ωr=0.85ω0)，橫搖角位移隨時(shí)間的變化曲線如圖7所示。建立絕對(duì)坐標(biāo)系，位于初始時(shí)刻箱體的下壁面中點(diǎn)，坐標(biāo)系不隨箱體振動(dòng)而運(yùn)動(dòng)。二維箱體初始水深h=0.093 m，內(nèi)部區(qū)域采用矩形布點(diǎn)，粒子間距Δx=Δy=0.002 m，CFL常數(shù)設(shè)置為0.2。

圖7 橫搖角位移曲線Fig.7 Time history of roll angle

為了驗(yàn)證劇烈晃蕩工況下模型的精度，在距離箱體底部0.093 m處的左側(cè)壁面設(shè)置一個(gè)壓力監(jiān)測(cè)點(diǎn)，坐標(biāo)值為B(x=-0.45 m,y=0.009 3 m)，用來(lái)檢測(cè)該點(diǎn)水壓隨時(shí)間的變化。圖8為采用兩種邊界粒子法計(jì)算出的B點(diǎn)處壓強(qiáng)隨時(shí)間變化的曲線以及與文獻(xiàn)[32]試驗(yàn)值的比較。由圖中可見(jiàn)兩種方法計(jì)算出的壓力峰值的時(shí)刻與試驗(yàn)值是一致的，三條壓力曲線的趨勢(shì)變化也基本一致，算例表明對(duì)于水壓值的計(jì)算，兩種方法均具有良好的計(jì)算精度。

圖8 點(diǎn)B處壓強(qiáng)比較Fig.8 Pressure comparison at point B

圖9為兩種邊界粒子法自由液面高度變化模擬結(jié)果與實(shí)驗(yàn)值的比較，可見(jiàn)兩種方法均能捕捉到艙壁抨擊、液體飛濺(t=3.45 s、t=4.35 s)以及自由液面翻卷(t=3.95 s、t=4.75 s)現(xiàn)象，適用于劇烈液體晃蕩現(xiàn)象數(shù)值模擬。

圖9 邊界粒子運(yùn)動(dòng)(圖左)、邊界粒子固定不動(dòng)(圖中)計(jì)算結(jié)果與試驗(yàn)值(圖右)比較Fig.9 Comparison of the numerical results of moving boundary particles(left)and fixed particles(middle) with experiments(right)

圖11為相同時(shí)間步長(zhǎng)不同粒子間距下計(jì)算結(jié)果比較，可見(jiàn)隨粒子間距的變小兩種方法計(jì)算值與試驗(yàn)值均誤差變小，邊界粒子固定法在粒子間距變化時(shí)，計(jì)算結(jié)果變化較小，相對(duì)于邊界粒子運(yùn)動(dòng)法，邊界粒子固定法具有更好的收斂性。

圖10 邊界粒子運(yùn)動(dòng)(圖左)、邊界粒子固定不動(dòng)(圖右)不同時(shí)間步長(zhǎng)下計(jì)算結(jié)果比較Fig.10 Comparison of the numerical results of moving boundary particles(left)and fixed particles(right) by different time steps

圖11 邊界粒子運(yùn)動(dòng)(圖左)、邊界粒子固定不動(dòng)(圖右)不同粒子間距下計(jì)算結(jié)果比較Fig.11 Comparison of the numerical results of moving boundary particles(left)and fixed particles(right)by different particle spacings

4 結(jié)論

以養(yǎng)殖工船養(yǎng)殖水艙晃蕩現(xiàn)象為研究對(duì)象，運(yùn)用SPH法分別對(duì)橫搖激勵(lì)下輕微液體晃蕩現(xiàn)象和劇烈液體晃蕩現(xiàn)象展開(kāi)模擬研究，并引入動(dòng)態(tài)邊界粒子法解決SPH法固壁邊界上或邊界附近的粒子截?cái)嗾`差影響精度問(wèn)題，實(shí)現(xiàn)了對(duì)SPH法的改進(jìn)。研究發(fā)現(xiàn)，邊界粒子運(yùn)動(dòng)法、邊界粒子固定不動(dòng)法均能有效地進(jìn)行液艙晃蕩現(xiàn)象數(shù)值模擬；相對(duì)于邊界粒子運(yùn)動(dòng)法，邊界粒子固定不動(dòng)法具有更好的算法穩(wěn)定性。本研究方法適用于液體晃蕩現(xiàn)象數(shù)值模擬，可以為養(yǎng)殖工船養(yǎng)殖水艙晃蕩模擬提供數(shù)值方法參考。

□