基于粒子法的規則波中浮體運動仿真

宋學敏 王緒明 劉維勤

(高性能艦船技術教育部重點實驗室1) 武漢 430063) (武漢理工大學船海與能源動力工程學院2) 武漢 430063)(武漢理工大學國家水運安全工程技術研究中心3) 武漢 430063)

0 引 言

無網格法基本理念是將連續系統以一系列運動粒子來近似，以拉格朗日技術描述空間運動粒子(遵循流體控制方程)，以運動粒子表征物理邊界且無需網格重構.由此避免了高海況下自由面大變形中網格畸變問題，適用于捕捉強非線性波浪兩相流問題(如破波、波浪翻卷等).與網格類方法相比，以拉格朗日法為核心技術的無網格法，計算中無需涉及網格結構，避免了網格生成、網格重構及網格質量所面臨的各種問題.

無網格法包括光滑粒子水動力學SPH方法和移動粒子半隱式法MPS方法.目前多數采用顯式弱可壓縮SPH方法[1-3]，包括采用多相流不可壓縮ISPH算法[4-6].MPS方法適用于研究自由液面大變形問題如自由液面流體流動問題、破波問題、波浪翻卷.越塚誠一等[7-8]開發了MPS方法并模擬了伴有流體碎片的不可壓縮粘性流體流動，柴田和也等[9-11]基于MPS方法研究了甲板上浪沖擊的影響，越塚誠一等[12]將MPS方法運用到流固耦合問題的研究中.求解多相流問題的關鍵技術涉及：交界面粒子的邊界捕捉技術；有效實施交界面粒子速度與壓力邊界條件；局部粒子堆積和分離的移動處理方法；多相壓力求解的解耦技術.傳統的MPS方法存在搜索粒子耗時、復雜邊界處理困難，易產生非物理壓力數值振蕩(如負壓顆粒間空隙、粒子分布不均一)，以及自由液面判斷準則方式(使某些純液體粒子誤判為自由液面粒子).經過近20年的發展，國內外學者從核函數、壓力梯度模型、Laplacian模型、壓力Poisson方程中的源項及自由面粒子的判斷精度等方面對MPS方法的計算穩定性、精度和效率都進行了系列的改進.Khayyer等[13]對壓力梯度的計算進行了改進，提出了一種動量守恒形式(即CMPS)，數值結果表明能在一定程度上緩解壓力的振蕩現象.隨后，Khayyer等[14]又提出修改壓力Poisson方程源項，提出了高階源項法(MPS-HS)，數值結果表明結合CMPS和MPS-HS方法，即CMPS-HS方法，能夠獲得較好的效果.萬德成等[15-17]基于改進MPS方法、動態負載平衡和MPI(message passing interface)庫自主開發了三維并行MPS方法求解器，并進一步將其應用于模擬自由面的流動問題.這些研究工作很大程度地推進了移動粒子半隱式法的發展，但從計算結果上看，這些改進的手段針對不同問題效果不一，還需要更多的應用來驗證.有些改進方法實施起來較為繁瑣，也需要提出解決問題的新思路.

文中采用改進的MPS方法進行線性波、Stokes二階波和浮體運動的數值模擬與分析，利用移動粒子半隱式法進行三維的線性波和二階Stokes波數值波浪水池的研究，基于線性波浪理論生成線性波，并計算結果和理論值分析對比，以及二階Stokes波的計算和對比驗證，驗證了MPS方法應用于波浪計算的可行性和可靠性.運用胡克定律模擬浮體的系泊系統，并進行規則波中浮體運動的數值仿真，分析浮式臺的運動響應.

1 移動粒子半隱式法

以半隱式算法為核心元素的MPS方法，通過Poisson方程隱式求解壓力，確保流體的不可壓縮性，顯著提升計算穩定性.圖1為傳統的MPS法流程框架示意圖.

圖1 傳統的MPS法流程框架示意圖

1.1 控制方程

在移動粒子半隱式法中，對于不可壓縮流體而言，其控制方程是質量和動量守恒方程，為

(1)

(2)

式中：ρ為流體的密度；u為流體的速度；p為壓力；μ為運動黏性系數；g為重力加速度.在移動粒子半隱式法中粒子的位置、速度等都是基于拉格朗日描述，粒子法的優勢是不用計算對流項.網格法用空間固定計算點來進行計算，速度的拉格朗日微分必須計算式(3)的右邊兩項，第二項是對流項.

(3)

1.2 粒子相互作用模型

在移動粒子半隱式法中，采用梯度模型和拉普拉斯模型進行控制方程的模擬，為

(4)

(5)

式中：d為空間的維度，二維空間時d=2，三維空間時d=3；n0為標準的粒子數密度；ri為第i個粒子的位置向量；φi為第i個粒子的任意的量；w(r)為由粒子間距離r和影響半徑re表示的核函數，此處采用經典的核函數模型考慮粒子間的相互影響見式(6).φi為第i個粒子的最小值，見式(7).

(6)

(7)

φvirtual=0

(8)

在移動粒子半隱式法中，式(5)為由越塚誠一等給出的拉普拉斯模型.其中引入了λ使得數值結果與擴散方程的解析解一致，λ為

(9)

1.3 改進的技術

MPS方法采用了半隱式算法求解N-S方程.N-S方程求解過程中，第一階段，采用顯示解法求解黏性項和重力項；第二階段，采用隱式解法求解泊松方程.

(10)

式中：pvirtual為大氣壓力.如果pi,pj≥0且pvirtual=0 Pa,式(10)就和式(11)等價

(11)

2 規則波的計算

2.1 流入流出算法

為了生成波浪，在水池的兩端使用了流入流出算法作為邊界條件，見圖2.由圖2可知，左右兩邊的兩列假想單元格分別為流入和流出邊界.外層的A層假想單元格用作流入單位，粒子會被生成，內側的B層假想單元格稱為輔助單元格用于計算所生成粒子的位置.等邊長的立方體表示假想的單元格，立方體的邊長等于粒子間的初始距離.流入邊界上，如果外側的A層假想單元格沒有粒子，粒子就會被生成在自由液面以下的空的單元格內.新生成的粒子被放在距離鄰近外側假想單元格高度和寬度方向上一定距離的地方.在流動方向的坐標中，這個距離被固定為初始的鄰近粒子間距.流出邊界上，一旦粒子流出了邊界的假想單元格，它們就會被刪除.有了流入和流出邊界，可以定義任意尺寸的數值波浪水池.

圖2 數值波浪水池示意圖

在水池中沿著Y軸等間距布置虛擬的圓柱體，測量進入虛擬圓柱體的流體粒子的最大Z軸坐標值，由此，測量得到生成的波浪時程曲線.

2.2 邊界上的壓力和速度

基于線性波浪理論的理論解析方法，根據伯努利定理的計算值設定邊界粒子的壓力，可以得到線性波的數值水池計算.同樣的，基于二階Stokes波浪理論的理論解析方法，根據伯努利定理的計算值設定邊界粒子的壓力，可以在數值水池中生成二階Stokes波.

方程(12)～(16)表達了二階Stokes理論.邊界流體粒子的X軸和Z軸速度分量由二階Stokes波浪理論的解析解給出.

(2+cosh 2kh)cos2(kx-ωt)

(12)

(13)

(14)

(15)

ω2=gktanhkh

(16)

2.3 數值水池的驗證

將數值水池中計算所得的結果與線性和二階Stokes波浪理論解析解的結果進行對比分析，從而驗證基于移動粒子半隱式法的數值波浪水池的可行性.計算區域的尺寸是6.72 m×0.6 m×1.00 m.驗證過程的波浪參數分別為：波長6.72 m、波高0.30 m和波浪周期2.075 s.一階和二階Stokes波浪的參數和流體速度的解析解作為初始值和邊界值根據波浪理論公式計算得到.在公式的計算中，水深h=8.00 m.實際上，為了節約計算時間，將計算區域可見的水池深度設為1.00 m，然而真實的波浪參數中水深為8.00 m.數值水池底部的粒子的X軸和Z軸的速度根據式(14)～(15)計算得出.

在此驗證計算中，流體的密度ρ為1 000 kg/m3.初始的粒子間距為0.03 m.初始時刻的粒子總數為283，552個.進行了20.75 s共計10個波浪周期的波浪傳播過程的計算.

2.4 數值結果的分析

通過與二階Stokes波的理論解進行對比分析，驗證了基于移動粒子半隱式法的波浪仿真計算的可行性.圖3為一個波浪周期內0.2 s時間間隔波浪粒子圖，波浪從右向左傳播，整個仿真過程中完整保持了波浪的振幅.圖4為MPS方法的線性波的波高時程曲線的計算結果和線性波的理論值的時程曲線的對比.粗曲線代表仿真計算結果，虛線代表線性波浪理論的理論值，由圖4可知：MPS方法的計算結果和線性波浪理論的解析值相當接近，仿真值和理論值最大的相差值是0.025 m相當于8.3%的波幅，這是工程應用中能夠接受的誤差范圍.圖5為MPS方法的二階Stokes波的波高時程曲線的計算結果和二階Stokes波的理論值的時程曲線的對比.同樣地，曲線代表仿真計算結果，虛線代表二階Stokes波浪理論的理論值，仿真值和理論值之間的差值非常小，仿真值和理論值最大的相差值是0.031 7 m相當于10.6%的波幅.通過上述結果的分析發現，移動粒子半隱式法不僅僅可以用于線性波浪的計算，同時也能用于二階Stokes波的計算.

圖3 各時刻波浪速度云圖

圖4 MPS方法的線性波的仿真結果和理論值的時間曲線

圖5 MPS方法的二階Stokes波的仿真結果波高和理論值的時間曲線

3 系泊的浮體運動仿真

3.1 系泊模型

假設系泊模型由極細的彈性繩索組成，可以忽略繩索與流體之間的相互作用.因此，可以用胡克定律來簡化計算系泊模型力.系泊模型的長度可以通過剛體上的固定點與水底的距離來計算.根據實際系泊模型，當長度小于初始長度時，系泊力Fm為零.相反，系泊力Fm為

(17)

根據式(18)～式(19)，計算剛體上系泊點a與水底系泊點b之間的不動點長度，(Xa,Ya,Za),(Xb,Yb,Zb)分別為“a”“b”點的x，y，z方向上的位置.系泊模型示意圖見圖6.初始長度是一個常數L0.

圖6 系泊模型

ΔL=L-L0

(18)

(19)

3.2 浮式平臺仿真

為了驗證本研究中數值方法的可靠性，進行了文獻[17]中水池試驗工況的數值模擬，具體模型參數、工況波浪參數和文獻中一致.運用三維設計軟件建立了浮式平臺的三維模型，利用三維軟件的數據生成等粒子間距的輸入文件，粒子間距等于0.01 m.數值模擬中采用與水池試驗相同的浮式平臺主要特征參數和相同的波浪參數.水池試驗情況見圖7.

圖7 水池試驗示意圖

根據文獻[17]中的試驗條件建立了規則波條件下的三維浮式平臺模型，并對該模型進行了數值模擬.在數值計算中：①在三維模型中考慮系泊模型；②建立了三維浮式平臺模型；③實現了浮式平臺6個自由度的運動的計算.通過對浮式平臺運動幅值的簡單平均得到MPS法的計算結果，表1中比較了浮式平臺的縱蕩、垂蕩、縱搖運動的振幅.結果表明，數值計算結果和試驗結果誤差在15%以內滿足工程應用的要求.

表1 數值計算結果和試驗結果對比

4 結 束 語

文中在傳統的MPS方法的基礎上加入了改進的壓力梯度模型，依托該改進的程序分別進行了線性波和二階Stokes波的數值波浪水池研究和模擬，通過MPS方法的數值計算結果和理論值結果對比分析，發現MPS方法可以生成任意波浪參數的線性波和二階Stokes波，并且計算結果滿足工程應用的精度要求.然后，基于改進的MPS方法的數值造波算法，引入了胡克定律模擬浮體的系泊系統，開展了規則波中浮體運動仿真.本研究一方面驗證了MPS方法應用于數值造波計算的有效性和可靠性，為研究航行中船舶砰擊、甲板上浪、船舶水彈性等強非線性流體計算提供了可能性.另一方面將該數值方法應用于波浪中系泊浮體運動的模擬，為MPS方法在船舶工程領域更廣泛的應用提供新思路.