基于改進(jìn)MPS方法的三維時(shí)域水彈性計(jì)算模型

孫 哲，尹恒輝，張桂勇，3，4，宗 智，3，4

（1.大連理工大學(xué)船舶工程學(xué)院 遼寧省深海浮動(dòng)結(jié)構(gòu)工程實(shí)驗(yàn)室，遼寧 大連116024；2.大連中遠(yuǎn)海運(yùn)重工有限公司，遼寧 大連116113；3.大連理工大學(xué) 工業(yè)裝備與結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連116024；4.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240）

0 引 言

與傳統(tǒng)的基于網(wǎng)格的CFD方法相比，粒子類方法以其拉格朗日形式的時(shí)間步進(jìn)格式和無網(wǎng)格的空間離散方法，為處理大自由液面變形問題提供了一個(gè)更加有效的計(jì)算模型。MPS（Moving Particle Semi-implicit）方法[1]和SPH（Smoothed Particle Hydrodynamics）方法[2]都是主流的粒子類方法。其中MPS方法最先基于投影法（Projection Method）[3]提出所謂半隱式求解格式，即通過求解壓力Poisson方程來實(shí)現(xiàn)壓力和速度的解耦，而非像在傳統(tǒng)弱可壓SPH（Weakly Compressible SPH）中通過狀態(tài)方程來基于密度變化求解壓力。這種半隱式求解方式可以使得時(shí)間步長(zhǎng)取得更大，并且能夠更好地保證不可壓縮性。當(dāng)然，原始MPS方法也和SPH一樣存在著壓力計(jì)算有較大非物理性震蕩問題，因此許多研究者對(duì)MPS進(jìn)行了改進(jìn)[4-8]，取得了較好的效果。

另一方面，由于現(xiàn)代船舶的大型化、高運(yùn)營(yíng)速度等發(fā)展趨勢(shì)，其固有頻率越來越接近于波浪作用的遭遇頻率，且受到波浪砰擊、甲板上浪等強(qiáng)非線性瞬態(tài)載荷的作用更加頻繁，會(huì)發(fā)生所謂砰擊顫振（Whipping）和波擊振蕩（Springing）等現(xiàn)象，嚴(yán)重地影響船舶的總體和疲勞強(qiáng)度[9]。這就要求在對(duì)結(jié)構(gòu)強(qiáng)度分析時(shí)采用水彈性方法，尤其是強(qiáng)非線性的時(shí)域水彈性計(jì)算是較為必要的。傳統(tǒng)的基于勢(shì)流的水彈性方法已經(jīng)形成了節(jié)點(diǎn)多、分支系的理論體系[10]。然而，在計(jì)及上述設(shè)計(jì)大自由面變形的強(qiáng)非線性時(shí)域水彈性計(jì)算方面，勢(shì)流理論還存在缺陷。事實(shí)上正如Hirdaris[11]在其綜述文章中所說的那樣，部分或全部使用CFD方法進(jìn)行水彈性計(jì)算將是終極目標(biāo)。本文將針對(duì)基于MPS方法的時(shí)域水彈性模型進(jìn)行探討。

本文首先闡述改進(jìn)MPS模型及其三維拓展，主要是三維情況下的自由面粒子搜索算法；其次將簡(jiǎn)要介紹基于拉格朗日方程的耦合剛體和彈性模態(tài)的結(jié)構(gòu)動(dòng)力響應(yīng)模型；最后，通過對(duì)三維潰壩和三維浮梁與潰壩問題的相互作用問題的模擬來驗(yàn)證模型的有效性。

1 改進(jìn)的MPS方法

本節(jié)將對(duì)作者所發(fā)展的改進(jìn)MPS方法進(jìn)行簡(jiǎn)要闡述，詳細(xì)算法流程參見文獻(xiàn)[6，12]。其中將重點(diǎn)介紹本文對(duì)該模型的三維拓展，即自由面粒子搜索算法。

1.1 半隱式求解格式

由于所研究問題粘性影響很小，因此流體部分控制方程為無粘不可壓流體NS方程和連續(xù)性方程，即：

式中：u，p和ρ分別為流體速度、壓強(qiáng)和密度；g為重力加速度矢量形式。

速度和壓強(qiáng)采用如下的投影法（Projection Method）來解耦：

首先，不計(jì)壓強(qiáng)，計(jì)算一個(gè)中間速度和對(duì)應(yīng)的中間位置，并將流體質(zhì)點(diǎn)步進(jìn)到該位置，即：

其次，在改進(jìn)MPS方法中，根據(jù)連續(xù)性方程構(gòu)造如下形式的壓力Poisson方程：

方程中n0和nn稱為“粒子密度”，其定義由（11）式給出，下標(biāo)0和n分別表示初始和上一時(shí)刻對(duì)應(yīng)的狀態(tài)。（3）式中系數(shù)α起到提高壓強(qiáng)計(jì)算穩(wěn)定性和準(zhǔn)確度的作用，其定義如下：

從上式可知，其計(jì)算不含人為的經(jīng)驗(yàn)假設(shè)。最終，粒子速度和位置可由計(jì)算得到的下一時(shí)刻的壓強(qiáng)來更新如下：

為了進(jìn)一步提高算法的穩(wěn)定性和精度，在改進(jìn)MPS方法中，粒子位置將根據(jù)分布情況被重新排布，具體細(xì)節(jié)可參見文獻(xiàn)[6]。

1.2 算子離散方法

在MPS方法中，梯度、散度和拉普拉斯算子是通過一種加權(quán)平均的形式離散，其具體格式如下：

式中：φ和φ分別表示任意標(biāo)量和矢量；d表示空間維度（即2或3維）；M表示在以中心粒子為圓心的圓形支撐域內(nèi)的粒子數(shù)目；對(duì)于梯度和散度算子，圓半徑re=2.1r0，對(duì)于拉普拉斯算子re=4.0r0，其中r0為粒子初始間距；w（rij）是權(quán)函數(shù)；λ是與權(quán)函數(shù)相關(guān)的參數(shù)。它們的定義如下：

前述的粒子密度定義如下：

由于粒子位置處于時(shí)刻變動(dòng)過程中，因此需要每個(gè)時(shí)間步搜索更新各個(gè)粒子的支持域內(nèi)粒子數(shù)目，在改進(jìn)MPS法中進(jìn)一步發(fā)展了一種基于背景規(guī)則網(wǎng)格（Cell-link）的搜索方式，所需搜索面積比原方法減少6.5/9。

1.3 壓力Poisson方程邊界條件

為求解壓力Poisson方程（即（3）式），需施加兩類邊界條件，即固體邊界條件和自由面邊界條件。

（a）固體邊界條件

在改進(jìn)MPS方法中，固體邊界上的粒子將被施加如下的Neumann邊界條件：

（b）自由面邊界條件

求解Poisson方程（即（3）式）需在自由液面處施加壓力為零的邊界條件，因此需搜索并確定自由液面粒子的位置。在作者所發(fā)展的二維改進(jìn)MPS算法中，自由液面粒子搜索采用對(duì)Koh等人[5]提出模型的改進(jìn)形式，即對(duì)每一個(gè)粒子分配一個(gè)圓環(huán)，如果圓環(huán)上均布的360個(gè)點(diǎn)被周圍粒子的圓環(huán)覆蓋，則認(rèn)為該粒子為內(nèi)部粒子，否則為自由液面粒子。本文中進(jìn)一步拓展該模型至三維情形，自由液面粒子判斷分為兩部分：

圖1 三維自由液面粒子判斷示意圖Fig.1 Sketch of the 3D free surface identification

首先，按照每個(gè)粒子圓形支撐域內(nèi)的粒子數(shù)目Np進(jìn)行第一輪判斷：

式中：Np0為均勻分布的情況下，典型內(nèi)部粒子的支撐域內(nèi)鄰近粒子數(shù)目，取為32；β3d為調(diào)節(jié)系數(shù)，本文中取為0.96。當(dāng)（13）式滿足時(shí)，該粒子被初步判斷為自由面粒子，否則為內(nèi)部粒子。

其次，為進(jìn)一步修正上一步中被誤判的粒子，為上一步中判斷出的自由液面粒子及其周圍粒子分別分配一個(gè)球面，并根據(jù)周圍粒子分布疏密程度，用權(quán)函數(shù)加權(quán)平均形式確定球面的法向方向（如圖1中的nz方向）。之后，以與法線方向夾角為θ的半徑，在球面上截取一個(gè)如圖1中紅色區(qū)域所示的圓形三維曲面，并在其上均布若干節(jié)點(diǎn)，如果這些節(jié)點(diǎn)全部被周圍粒子的球面覆蓋，則認(rèn)為該中心粒子為內(nèi)部粒子，否則為自由面粒子。

這樣的搜索方式的好處是結(jié)合了粒子數(shù)判別法的高效率和基于幾何分布判別法的準(zhǔn)確性，下面將給出的例子證明了所提出模型的有效性。

2 時(shí)域耦合剛體和彈性模態(tài)計(jì)算模型

對(duì)于在波浪中運(yùn)動(dòng)船舶來說，最為重要的變形就是縱向彎曲變形，其核心是一個(gè)二維變形過程（縱向剖面內(nèi)）。作者基于拉格朗日方程，提出了一個(gè)耦合剛性和彈性模態(tài)的結(jié)構(gòu)響應(yīng)二維計(jì)算模型[13]，可用于在時(shí)域高效地計(jì)算大剛體位移疊加相對(duì)較小彈性變形的典型船體運(yùn)動(dòng)情形。下面簡(jiǎn)要描述該模型。

2.1 運(yùn)動(dòng)描述

圖2 二維彈性浮梁模型示意圖Fig.2 Sketch of the 2D floating beam

根據(jù)上述的坐標(biāo)系統(tǒng)，梁上任一點(diǎn)的位置X可由下式計(jì)算：

式中：R為用于貼體坐標(biāo)系和全局坐標(biāo)系之間轉(zhuǎn)換的旋轉(zhuǎn)矩陣；ξ為變形后的點(diǎn)在貼體坐標(biāo)系中的表達(dá)式。R和ξ分別由（15）式和（16）式定義：

式中w為梁未變形時(shí)任意一點(diǎn)對(duì)應(yīng)的貼體坐標(biāo)系中O-w軸方向的坐標(biāo)值。梁上任一點(diǎn)的速度和加速度可由下式求得：

2.2 模態(tài)疊加

應(yīng)用模態(tài)疊加理論，將撓度部分表示成如下形式：

式中Φ和q為各階模態(tài)函數(shù)向量以及對(duì)應(yīng)的廣義坐標(biāo)向量，它們的定義如下：

通過這種模態(tài)展開的形式實(shí)現(xiàn)了時(shí)間和空間的變量分離。對(duì)（19）式求一階和二階時(shí)間導(dǎo)數(shù)得到：

模態(tài)函數(shù)之間滿足如下的正交關(guān)系：

式中：ωk表示梁的第k階振動(dòng)固有圓頻率；E，ρl和J分別是材料楊氏模量、線密度和截面二階慣性矩；I是n階單位矩陣；n是所取得彈性模態(tài)數(shù)量。式中的線積分是沿著梁的中心線進(jìn)行的。在本文的研究中，應(yīng)用的是自由-自由梁的模態(tài)函數(shù)：

式中：M為梁的總質(zhì)量；a=L/2為梁的長(zhǎng)度一半值；μk為下式的正實(shí)數(shù)特征值：

前三階取值為：μ1=2.365 0，μ2=3.926 6，μ3=5.497 8。

2.3 基于拉格朗日方程的模型

根據(jù)上述運(yùn)動(dòng)描述模型，梁的變形可由以下廣義位置變量D和對(duì)應(yīng)的廣義力Q描述：

式中廣義力變量QXcR，Qθ和Qqi是對(duì)應(yīng)的剛體（即Xc和θ）和彈性（即q）部分的非保守力，可由（30）～（32）式計(jì)算：

式中：n表示結(jié)構(gòu)表面任意一點(diǎn)的外法線方向向量（指向結(jié)構(gòu)內(nèi)部）；ew表示貼體坐標(biāo)系中O-w軸方向的單位向量。根據(jù)拉格朗日方程，結(jié)構(gòu)的動(dòng)力響應(yīng)可表示為：

式中T和V分別為結(jié)構(gòu)的動(dòng)能和勢(shì)能：

將結(jié)構(gòu)動(dòng)能（（36）式）和勢(shì)能（（37）式）代入拉格朗日方程（（33）～（35）式）中，并取前三階彈性模態(tài)，化簡(jiǎn)得如下浮梁在波浪中結(jié)構(gòu)動(dòng)態(tài)響應(yīng)的控制方程：

式中：Iin為梁圍繞其質(zhì)心的二階轉(zhuǎn)動(dòng)慣性矩（非截面慣性矩）；ψ0i和ψ1i（i=1,2,…n）定義如下：

（38）～（43）式所組成的非線性二階微分方程組中考慮了剛體和彈性模態(tài)的相互影響，即交叉項(xiàng)。該方程組將通過Newmark方法[14]降階為一階線性方程組后由Newton-Raphson方法來進(jìn)行求解。流固耦合部分采用基于Aitken松弛模型的迭代方式求解，具體細(xì)節(jié)可參見文獻(xiàn)[13]。

3 計(jì)算結(jié)果及分析

3.1 改進(jìn)MPS模型驗(yàn)證計(jì)算—3D潰壩問題

為檢驗(yàn)三維改進(jìn)MPS模型的有效性，本文對(duì)圖3所示的三維潰壩與方形平臺(tái)相互作用模型進(jìn)行了模擬，這可以認(rèn)為是甲板上浪問題的近似模型。若干水位（H1-H4）和壓力（p1-p8）傳感器被布置在了如圖3中所示的水池和和池中的方形平臺(tái)上。模擬中采用的初始粒子間距為0.015 m，共308 460個(gè)粒子，其中200 244個(gè)流體粒子。時(shí)間步長(zhǎng)根據(jù)CFL條件動(dòng)態(tài)調(diào)整，初始值為0.002 s。

典型位置處的水位和壓力監(jiān)測(cè)曲線與實(shí)驗(yàn)和商業(yè)軟件ComFLOW計(jì)算結(jié)果[15]的對(duì)比如圖4所示，其中圖4（a）～（b）表示圖3中H2和H4位置處的水位時(shí)程曲線，圖4（c）～（d）表示圖3中p3和p7位置處的壓力時(shí)程曲線。從圖中結(jié)果可以看出，改進(jìn)的三維MPS計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好，且并未出現(xiàn)ComFLOW計(jì)算結(jié)果中的奇異值現(xiàn)象（約t=1.4 s處）。

圖3 三維潰壩模型示意圖Fig.3 Sketch of the 3D dam-break model

圖5給出了自由面形態(tài)和實(shí)驗(yàn)以及商業(yè)軟件ComFLOW計(jì)算結(jié)果的對(duì)比，左側(cè)MPS結(jié)果中僅用紅色標(biāo)記出自由面粒子的位置，可以看出在劇烈的變形過程中并未出現(xiàn)粒子誤判現(xiàn)象，證明了所發(fā)展的三維自由面粒檢測(cè)方法是有效的。兩個(gè)典型的自由面形態(tài)與實(shí)驗(yàn)及ComFLOW結(jié)果相比吻合良好。

3.2 3D潰壩沖擊彈性浮梁計(jì)算

圖4 三維潰壩典型位置處水位和壓力時(shí)程曲線（上：水位；下：壓力）Fig.4 Time history of water level and pressure(Up：Water level;Bottom：Pressure)

本節(jié)將對(duì)一個(gè)二維變形的浮梁在三維潰壩沖擊作用下的整體動(dòng)力響應(yīng)進(jìn)行計(jì)算，來驗(yàn)證本文給出的流固耦合模型對(duì)此類強(qiáng)非線性時(shí)域水彈性問題的有效性。計(jì)算模型布置如圖6所示。計(jì)算中所用的初始粒子間距為0.02 m，共194 360個(gè)粒子，其中流體粒子數(shù)為126 288。時(shí)間步長(zhǎng)同樣以CFL為限制條件動(dòng)態(tài)控制，初始值是0.002 s。為驗(yàn)證模型在處理較大的彈性變形的有效性，方梁的剛度取得很小，第二節(jié)中的結(jié)構(gòu)計(jì)算模型的參數(shù)EJ=4.5 N/m2，線密度ρl=36 kg/m，對(duì)應(yīng)的梁的前三階固有圓頻率為ω1=7.910 0 rad/s，ω2=21.804 6 rad/s，ω3=42.745 7 rad/s。

圖5 三維潰壩問題自由液面形態(tài)對(duì)比（上：t=0.4 s；下：t=0.56 s）Fig.5 Free surface profiles comparison(Up：t=0.4 s;Bottom：t=0.56 s)

圖6 三維潰壩沖擊浮梁模型示意圖Fig.6 Sketch of the 3D dam-break impacting floating beam model

圖7 典型時(shí)刻自由面形態(tài)及浮梁變形（（a）t=0.1 s；（b）t=0.418 s；（c）t=0.864 s；（d）t=1.282 s）Fig.7 Free surface profiles and beam deflections at typical time instants((a)t=0.1 s;(b)t=0.418 s;(c)t=0.864 s;(d)t=1.282 s)

從圖7中可以看出，在若干典型時(shí)刻的自由面粒子被正確地檢出，并且梁變形和流體的形態(tài)也是符合物理實(shí)際的。進(jìn)一步的系統(tǒng)化數(shù)值實(shí)驗(yàn)將在下一步的工作中開展，以全面驗(yàn)證模型的準(zhǔn)確性和可靠性。

4 結(jié) 論

本文首先將作者發(fā)展的二維改進(jìn)MPS方法進(jìn)行了三維拓展，主要是三維自由液面粒子的高效和準(zhǔn)確檢索模型。為驗(yàn)證模型的準(zhǔn)確性，對(duì)三維潰壩與平臺(tái)的沖擊問題進(jìn)行了模擬。數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果和其他數(shù)值模型吻合得較好，證明了所發(fā)展模型的有效性。在此基礎(chǔ)上，應(yīng)用作者所發(fā)展的耦合剛體和彈性模態(tài)的時(shí)域計(jì)算模型，對(duì)三維潰壩沖擊彈性浮梁?jiǎn)栴}進(jìn)行了模擬，得到的結(jié)果符合物理實(shí)際。下一步的工作將針對(duì)該水彈性模型進(jìn)行系統(tǒng)的數(shù)值驗(yàn)證，以證明其在處理此類涉及強(qiáng)非線性時(shí)域水彈性問題的有效性。