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

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

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

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

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

1 SPH法數(shù)值模型

1.1 控制方程及其離散形式

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

(1)

(2)

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

(3)

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

(4)

(5)

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

(6)

1.2 核函數(shù)

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

(7)

式中；對于二維問題αd=7/4πh2。

1.3 時間積分

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

在預估步有：

(8)

在校正步有：

(9)

(10)

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

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

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

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

2.2 動態(tài)邊界粒子運動

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

uw=u0+ω×(rw-rO)

(11)

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

在預估步，邊界粒子的位置為：

(12)

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

(13)

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

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

以物體壁面隨體坐標系為參考坐標系，保持壁面邊界不動，流場所受重力方向發(fā)生變化，對于二維問題，有如下公式：

(14)

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

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

(15)

式中：第一項為角加速度項，第二項為向心加速度項，第三項為科氏加速度項。

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

(16)

3 數(shù)值算例

3.1 輕微液體晃蕩算例分析

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

圖2 矩形箱體坐標系設置示意圖Fig.2 Schematic diagram of sloshing tank and the co-ordinate system

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

圖3 點A處壓強比較Fig.3 Pressure comparison at point A

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

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

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

圖5為兩種方法在模擬過程中幾個典型時刻壓力場的對比，左邊為邊界粒子運動的壓力場，右邊為邊界粒子不動的壓力場?？傮w來看兩種方法的模擬結(jié)果接近，均能很好地捕捉自由液面形狀。

圖5 計算結(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對箱體施加橫搖激勵使其繞箱體中心(軸x=0,y=0)做簡諧振動，根據(jù)箱體自振頻率表達式：

圖6 矩形箱體坐標系設置示意圖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)，橫搖角位移隨時間的變化曲線如圖7所示。建立絕對坐標系，位于初始時刻箱體的下壁面中點，坐標系不隨箱體振動而運動。二維箱體初始水深h=0.093 m，內(nèi)部區(qū)域采用矩形布點，粒子間距Δx=Δy=0.002 m，CFL常數(shù)設置為0.2。

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

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

圖8 點B處壓強比較Fig.8 Pressure comparison at point B

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

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

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

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

圖11 邊界粒子運動(圖左)、邊界粒子固定不動(圖右)不同粒子間距下計算結(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)象為研究對象，運用SPH法分別對橫搖激勵下輕微液體晃蕩現(xiàn)象和劇烈液體晃蕩現(xiàn)象展開模擬研究，并引入動態(tài)邊界粒子法解決SPH法固壁邊界上或邊界附近的粒子截斷誤差影響精度問題，實現(xiàn)了對SPH法的改進。研究發(fā)現(xiàn)，邊界粒子運動法、邊界粒子固定不動法均能有效地進行液艙晃蕩現(xiàn)象數(shù)值模擬；相對于邊界粒子運動法，邊界粒子固定不動法具有更好的算法穩(wěn)定性。本研究方法適用于液體晃蕩現(xiàn)象數(shù)值模擬，可以為養(yǎng)殖工船養(yǎng)殖水艙晃蕩模擬提供數(shù)值方法參考。

□