基于改進MPS方法的三維時域水彈性計算模型

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

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

0 引 言

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

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

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

1 改進的MPS方法

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

1.1 半隱式求解格式

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

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

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

首先，不計壓強，計算一個中間速度和對應的中間位置，并將流體質點步進到該位置，即：

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

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

從上式可知，其計算不含人為的經驗假設。最終，粒子速度和位置可由計算得到的下一時刻的壓強來更新如下：

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

1.2 算子離散方法

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

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

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

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

1.3 壓力Poisson方程邊界條件

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

（a）固體邊界條件

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

（b）自由面邊界條件

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

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

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

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

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

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

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

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

2.1 運動描述

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

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

式中：R為用于貼體坐標系和全局坐標系之間轉換的旋轉矩陣；ξ為變形后的點在貼體坐標系中的表達式。R和ξ分別由（15）式和（16）式定義：

式中w為梁未變形時任意一點對應的貼體坐標系中O-w軸方向的坐標值。梁上任一點的速度和加速度可由下式求得：

2.2 模態(tài)疊加

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

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

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

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

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

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

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

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

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

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

式中：n表示結構表面任意一點的外法線方向向量（指向結構內部）；ew表示貼體坐標系中O-w軸方向的單位向量。根據(jù)拉格朗日方程，結構的動力響應可表示為：

式中T和V分別為結構的動能和勢能：

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

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

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

3 計算結果及分析

3.1 改進MPS模型驗證計算—3D潰壩問題

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

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

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

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

3.2 3D潰壩沖擊彈性浮梁計算

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

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

圖5 三維潰壩問題自由液面形態(tài)對比（上：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 典型時刻自由面形態(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中可以看出，在若干典型時刻的自由面粒子被正確地檢出，并且梁變形和流體的形態(tài)也是符合物理實際的。進一步的系統(tǒng)化數(shù)值實驗將在下一步的工作中開展，以全面驗證模型的準確性和可靠性。

4 結 論

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