基于粒子法的規(guī)則波中浮體運(yùn)動(dòng)仿真

宋學(xué)敏 王緒明 劉維勤

(高性能艦船技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室1) 武漢 430063) (武漢理工大學(xué)船海與能源動(dòng)力工程學(xué)院2) 武漢 430063)(武漢理工大學(xué)國(guó)家水運(yùn)安全工程技術(shù)研究中心3) 武漢 430063)

0 引 言

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

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

文中采用改進(jìn)的MPS方法進(jìn)行線性波、Stokes二階波和浮體運(yùn)動(dòng)的數(shù)值模擬與分析，利用移動(dòng)粒子半隱式法進(jìn)行三維的線性波和二階Stokes波數(shù)值波浪水池的研究，基于線性波浪理論生成線性波，并計(jì)算結(jié)果和理論值分析對(duì)比，以及二階Stokes波的計(jì)算和對(duì)比驗(yàn)證，驗(yàn)證了MPS方法應(yīng)用于波浪計(jì)算的可行性和可靠性.運(yùn)用胡克定律模擬浮體的系泊系統(tǒng)，并進(jìn)行規(guī)則波中浮體運(yùn)動(dòng)的數(shù)值仿真，分析浮式臺(tái)的運(yùn)動(dòng)響應(yīng).

1 移動(dòng)粒子半隱式法

以半隱式算法為核心元素的MPS方法，通過(guò)Poisson方程隱式求解壓力，確保流體的不可壓縮性，顯著提升計(jì)算穩(wěn)定性.圖1為傳統(tǒng)的MPS法流程框架示意圖.

圖1 傳統(tǒng)的MPS法流程框架示意圖

1.1 控制方程

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

(1)

(2)

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

(3)

1.2 粒子相互作用模型

在移動(dòng)粒子半隱式法中，采用梯度模型和拉普拉斯模型進(jìn)行控制方程的模擬，為

(4)

(5)

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

(6)

(7)

φvirtual=0

(8)

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

(9)

1.3 改進(jìn)的技術(shù)

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

(10)

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

(11)

2 規(guī)則波的計(jì)算

2.1 流入流出算法

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

圖2 數(shù)值波浪水池示意圖

在水池中沿著Y軸等間距布置虛擬的圓柱體，測(cè)量進(jìn)入虛擬圓柱體的流體粒子的最大Z軸坐標(biāo)值，由此，測(cè)量得到生成的波浪時(shí)程曲線.

2.2 邊界上的壓力和速度

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

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

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

(12)

(13)

(14)

(15)

ω2=gktanhkh

(16)

2.3 數(shù)值水池的驗(yàn)證

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

在此驗(yàn)證計(jì)算中，流體的密度ρ為1 000 kg/m3.初始的粒子間距為0.03 m.初始時(shí)刻的粒子總數(shù)為283，552個(gè).進(jìn)行了20.75 s共計(jì)10個(gè)波浪周期的波浪傳播過(guò)程的計(jì)算.

2.4 數(shù)值結(jié)果的分析

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

圖3 各時(shí)刻波浪速度云圖

圖4 MPS方法的線性波的仿真結(jié)果和理論值的時(shí)間曲線

圖5 MPS方法的二階Stokes波的仿真結(jié)果波高和理論值的時(shí)間曲線

3 系泊的浮體運(yùn)動(dòng)仿真

3.1 系泊模型

假設(shè)系泊模型由極細(xì)的彈性繩索組成，可以忽略繩索與流體之間的相互作用.因此，可以用胡克定律來(lái)簡(jiǎn)化計(jì)算系泊模型力.系泊模型的長(zhǎng)度可以通過(guò)剛體上的固定點(diǎn)與水底的距離來(lái)計(jì)算.根據(jù)實(shí)際系泊模型，當(dāng)長(zhǎng)度小于初始長(zhǎng)度時(shí)，系泊力Fm為零.相反，系泊力Fm為

(17)

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

圖6 系泊模型

ΔL=L-L0

(18)

(19)

3.2 浮式平臺(tái)仿真

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

圖7 水池試驗(yàn)示意圖

根據(jù)文獻(xiàn)[17]中的試驗(yàn)條件建立了規(guī)則波條件下的三維浮式平臺(tái)模型，并對(duì)該模型進(jìn)行了數(shù)值模擬.在數(shù)值計(jì)算中：①在三維模型中考慮系泊模型；②建立了三維浮式平臺(tái)模型；③實(shí)現(xiàn)了浮式平臺(tái)6個(gè)自由度的運(yùn)動(dòng)的計(jì)算.通過(guò)對(duì)浮式平臺(tái)運(yùn)動(dòng)幅值的簡(jiǎn)單平均得到MPS法的計(jì)算結(jié)果，表1中比較了浮式平臺(tái)的縱蕩、垂蕩、縱搖運(yùn)動(dòng)的振幅.結(jié)果表明，數(shù)值計(jì)算結(jié)果和試驗(yàn)結(jié)果誤差在15%以?xún)?nèi)滿足工程應(yīng)用的要求.

表1 數(shù)值計(jì)算結(jié)果和試驗(yàn)結(jié)果對(duì)比

4 結(jié) 束 語(yǔ)

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