MPS-DEM耦合方法的多段系泊浮體波浪響應分析

孫一頡,席光,孫中國

(西安交通大學能源與動力工程學院,710049,西安)

海洋蘊含豐富的能量和資源,海洋資源和能源的開發(fā)離不開各種各樣的海上設施,其中海上浮體平臺作為一種常見的海上設施,在海洋工程中發(fā)揮著重要角色。為了保證大型浮體平臺運動的穩(wěn)定性和平臺上的人員安全,平臺結構常采用分段多模塊式連接。海浪中多模塊分段式耦合運動與單一整體自由運動存在明顯差異,了解分段式平臺的運動以及多模塊之間的動態(tài)響應關系對合理設計大型浮體平臺的結構十分重要。

海上浮體平臺與海浪的相互作用是一個復雜的強非線性問題,包含浮體結構受到的非線性水動力載荷和非線性變形。浮體平臺在海浪的沖擊下會發(fā)生變形,平臺的反復形變不僅可能會導致結構疲勞失效,還會引起海水與平臺相互作用力的變化,這類復雜問題通常采用無網(wǎng)格數(shù)值方法對浮體平臺的運動響應與結構受力進行預測。無網(wǎng)格法將計算區(qū)域離散為一系列的粒子,通過求解粒子的速度、位置、壓力及時間層推進來獲得流場的動態(tài)信息,與傳統(tǒng)網(wǎng)格方法相比,無網(wǎng)格法避免了網(wǎng)格重構和畸變等問題,在捕捉海浪破碎、甲板上浪等流體大變形現(xiàn)象和模擬固體變形斷裂等問題時有較大優(yōu)勢。

近年來,許多學者將無網(wǎng)格法或無網(wǎng)格法與其他方法的耦合算法應用于海上浮體結構物在海浪作用下的運動響應。Shibata等建立了三維剛性船體模型,并利用耦合移動粒子半隱式法(MPS)方法數(shù)值模擬了船舶在高海浪下的運動響應[1];Sun等采用可壓縮SPH方法模擬了剛體和波浪自由面的相互作用[2];Zhang等利用MPS-FEM方法,模擬了柔性浮體和波浪的相互作用[3];Sun等將MPS方法與模態(tài)疊加法耦合,計算了潰壩沖擊下浮體的變形與運動[4]。多數(shù)研究采用無網(wǎng)格法與傳統(tǒng)網(wǎng)格法相耦合的方法來計算海浪與浮體間的相互作用,使用傳統(tǒng)網(wǎng)格法離散的浮體結構區(qū)域仍存在網(wǎng)格,當浮體出現(xiàn)大變形或斷裂時網(wǎng)格易出現(xiàn)畸變,影響計算結果精度。

因此,本文耦合了MPS、離散單元法(DEM)兩種無網(wǎng)格方法,完全擺脫了數(shù)值計算對網(wǎng)格的依賴,充分發(fā)揮無網(wǎng)格法解決海浪與浮體非線性相互作用的優(yōu)勢,并利用MPS-DEM耦合算法,采用MPS方法模擬海浪運動,利用DEM方法離散分段浮體平臺。首先對二維活塞式造波模型進行了修正,獲得更加規(guī)則的波浪便于準確分析浮體響應;其次改進了消波模型并與解析解對比驗證了二維數(shù)值水槽造波的準確性;最后數(shù)值模擬了分段浮體平臺在規(guī)則波沖擊下的運動響應,分析了子平臺間的動態(tài)關系,討論了子平臺縱蕩、垂蕩和縱搖這3類運動響應。

1 數(shù)值方法

1.1 MPS方法簡介

MPS方法用于模擬液體的非定常大變形運動,液體粒子在拉格朗日框架下的控制方程為

(1)

(2)

式中:ρ為流體密度;μ為動力黏性系數(shù);u為流體速度;p為壓力;f為體積力。

粒子之間的相互作用強弱由以下核函數(shù)衡量

(3)

式中:w(r)為核函數(shù)值;r=|ri-rj|為任意兩個粒子i、j之間的距離;re為粒子相互作用半徑。梯度算子采用文獻[5]改進的模型,即

〈φ〉

(4)

壓力泊松方程計算公式改進[6]為

(5)

參數(shù)α取值范圍為(0,1),為了有效減小壓力計算誤差[7],本文取值為0.008。

本文液體自由表面的處理方式是將粒子數(shù)密度滿足n<βn0(β=0.97)條件的粒子判定為表面粒子,并將液體表面粒子壓力取值為0。

1.2 DEM方法簡介

DEM方法用于計算固體的運動、變形及相互作用,其中關于固體運動的流固耦合模擬,固體粒子運動控制方程[8]為

(6)

(7)

式中:m、I、u、ω分別為每個固體粒子的質(zhì)量、轉(zhuǎn)動慣量、速度和角速度;Fc、Tc分別為固體粒子間相互接觸產(chǎn)生的力和力矩;Fo、To分別為固體粒子受到的外力和外力矩;Fls、Tls為液體粒子對固體粒子的力和力矩。流體對固體的作用是通過壓力實現(xiàn)的,壓力作用于固體表面,方向沿固體表面法線方向指向固體內(nèi)部,相關計算公式為

Fls=-?SpdS=-?V

(8)

Tls=-?SrralpdS=-?VrralpdV=

(9)

式中:rral為粒子對于固體質(zhì)心的相對位置;l0為粒子直徑。

關于固體的變形模擬,即固體粒子間的相互作用與變形關系,離散單元法接觸模型如圖1所示。粒子間接觸力由剛性項和阻尼項構成,其法向分量和切向分量計算公式[9]為

Fcn=-knδn-dnun

(10)

(11)

式中:k、d、δ、u分別為剛度系數(shù)、阻尼系數(shù)、接觸位移和粒子速度;f為粒子間摩擦系數(shù);下標n表示法向;t表示切向。

圖1 DEM接觸模型

(a)粒子-粒子連接

(b)平行鍵連接圖2 DEM粒子連接模型

2 造波數(shù)值模型及改進

2.1 造波修正與精確造波

造波池造波有兩種方式:①直接設定造波板正弦運動[12];②由目標波形反演造波板運動規(guī)律[13]。第1種方式簡單易行,直接給定造波板的運動幅值與周期即可,但此方法造出波浪的實際參數(shù)與造波板運動規(guī)律并不完全一致,不適用于對波形參數(shù)要求較高的情況。斜坡式消波的活塞造波池幾何模型如圖3所示,粒子初始間距為0.025 m,粒子總數(shù)為16 095,水池左側(cè)為活塞式造波板,右側(cè)布置斜坡式消波段,防止反射波對波浪運動的影響。采用第1種方式進行造波,給定造波板運動

(12)

式中:l為造波板振幅;T為造波板運動周期。

圖3 斜坡式消波的活塞造波池幾何模型

直接造波法的自由面與解析解比較和反演方式造波的自由面與解析解比較如圖4、圖5所示。與振幅為0.1 m、周期為1.4 s的正弦波相比,雖然波浪周期與造波板運動周期相同,但是造波板振幅為0.1 m時,造出的實際波高更大且存在波動。圖5中采用反演方式造出的波形振幅和周期與解析解基本吻合。

圖4 直接造波法的自由面與解析解比較

圖5 反演方式造波的自由面與解析解比較

采用第2種方式造波時,活塞式造波板的水平方向速度運動規(guī)律為

(13)

(14)

ηp=Hcos(ωt)

(15)

式中:ω為波浪頻率;W為水力傳遞函數(shù);ηp為目標波形;k為波數(shù);d為水深;H為波高。

兩種造波方式下造波板的運動速度和位移比較如圖6所示。由圖6可知,在相同波形輸出條件要求下,反演式造波板的運動幅值和速度均小于直接造波。直接式造波方式雖然簡單,但是造出的波振幅被放大,而反演式造波通過目標波形條件來控制造波板運動,實現(xiàn)波浪參數(shù)實時反饋修正,適合應用于波浪條件已知、對波形參數(shù)要求較高的情況。

(a)造波板運動速度比較

(b)造波板橫向位移比較圖6 兩種造波方式造波板運動速度和位移比較

2.2 消波段改進

斜坡式消波段通過將液體橫向運動轉(zhuǎn)為豎直方向運動來逐漸降低波浪傳播速度,從而消除二次反射波,阻尼式消波的活塞造波池幾何模型如圖7所示。斜坡式消波效果較好,但在波長較長時,消波段所需長度較長,會大大增加計算量。本文采用阻尼式消波方式[13],在相同粒子尺寸時,粒子總數(shù)較斜坡式消波明顯降低,相對減少了12.1%,有效提高了計算效率。

圖7 阻尼式消波的活塞造波池幾何模型

消波段內(nèi)采用對粒子加速度的修正實現(xiàn)對波浪能的衰減,修正公式[13]為

(16)

阻尼式消波自由面高度與解析解比較如圖8所示,可知消波段改進后液面高度的監(jiān)測值和解析解吻合較好,說明消波形式的改變并未影響到池內(nèi)波形傳播,消波效果與斜坡式消波基本相同。

圖8 阻尼式消波自由面高度與解析解比較

3 算例與結果分析

本文基于MPS-DEM耦合方法,建立了雙向流固耦合模型,并通過數(shù)值計算與其他方法和實驗做了對比[8],驗證了方法的準確性。在此基礎上,數(shù)值模擬了分段浮體平臺在規(guī)則波沖擊下的運動,分段浮體平臺的首段和尾段分別系泊于造波池底部,不計系泊纜繩的質(zhì)量,如圖9所示。相鄰兩個子平臺間設立連接部件,平臺相對連接部件的運動有3個自由度,連接部件與子平臺內(nèi)設置連接粒子,連接粒子間作用包含彈性作用和阻尼作用,通過改變連接粒子間的彈性系數(shù)與阻尼系數(shù),達到控制平臺最大旋轉(zhuǎn)角度與最大平移的目的,平臺間約束作用示意圖如圖10所示。

圖9 分段浮體平臺幾何模型

(a)平臺間約束作用

(b)連接粒子間作用關系圖10 平臺間約束作用示意圖

二維情況下每段子平臺相對連接部件有縱蕩、垂蕩和縱搖3個自由度。本文針對含有2～5段子平臺的4種分段浮體平臺在波浪中的運動進行了數(shù)值模擬,并分析子平臺間的運動響應。通過分析不同系泊方式對浮體運動的影響[14],本文中的4組算例選取了定位效果更好的首尾雙纜繩系泊形式。分段浮體平臺算例參數(shù)如表1所示。

表1 分段浮體平臺算例參數(shù)

含有5節(jié)子平臺的浮體在規(guī)則波沖擊下運動的局部時序圖如圖11所示,深色部分為分段浮體平臺,淺色部分為海水。由圖11可知:t=2.0 s時規(guī)則波剛到達平臺前端,首段子平臺前仰帶動后部平臺;t=2.2 s時波浪沖擊首段子平臺,由于子平臺前端系泊,纜繩牽拉首段子平臺,平臺仰角減小浸入海水,少量海水躍上平臺;t為2.4～2.6 s時分段平臺主體受到波浪沖擊,分段結構呈彎曲狀態(tài),每個子平臺下部與海水緊密接觸,削弱了平臺主體的整體轉(zhuǎn)動;t為2.8～3.0 s時,波浪離開平臺,平臺恢復初始平展狀態(tài)。

(a)t=2.0 s (b)t=2.2 s

(c)t=2.4 s (d)t=2.6 s

(e)t=2.8 s (f)t=3.0 s圖11 5段浮體在規(guī)則波沖擊下運動的局部時序圖

(a)2段浮體

(b)3段浮體

(c)4段浮體

圖12為含有2～5節(jié)子平臺的4種浮體平臺的縱蕩響應,縱蕩指平臺沿最長延展方向的水平線性振蕩運動,即平臺沿x軸的振蕩。由于每個子平臺長度相同,因此相鄰子平臺間縱蕩曲線起始點(t=0 s)間隔相同;子平臺位置越靠近浮體尾部,縱蕩振幅越大,具體數(shù)據(jù)如表2所示,尾部的子平臺在波浪沖擊和連接部件約束下呈正弦運動,正弦運動幅值越大,說明尾部子平臺較前部子平臺需要連接部件更大的約束力以保證連接穩(wěn)定。在極端海浪沖擊下,大型分段浮體平臺需注意尾部子平臺水平方向的振蕩,防止子平臺間的相互碰撞和子平臺脫離分段浮體平臺主體。

(d)5段浮體圖12 4種浮體平臺的縱蕩響應

(a)2段浮體

(b)3段浮體

(c)4段浮體

(d)5段浮體圖13 4種浮體平臺的垂蕩響應

4種浮體平臺的垂蕩響應如圖13所示。由圖13可知:子平臺沿波浪傳播方向的位置先后,使得垂蕩曲線存在一定相位差;隨著波浪沿平臺尾部傳遞,能量逐段遞減,尾段子平臺垂蕩幅度最小;t=2.5 s時第2段子平臺垂蕩響應遠大于首段子平臺,浮體平臺在一個周期內(nèi)的平均運動響應如表2所示,這是由于首段子平臺前端系泊,當海浪浪高較大時系泊纜繩會對子平臺產(chǎn)生位移約束,避免子平臺運動幅值過大而導致分段浮體平臺主體振蕩不穩(wěn)。

(a)2段浮體

(b)3段浮體

(c)4段浮體

4種浮體平臺與水平線的夾角變化,反映了各個子平臺的縱搖響應,4種浮體平臺的縱搖響應如圖14所示,縱搖是平臺繞橫軸的回轉(zhuǎn)振蕩運動,文中定義平臺下俯時,與水平線的夾角為正。由圖14可知,波浪剛沖擊到首段子平臺時,平臺前端上仰,縱搖角度急劇下降由正轉(zhuǎn)負,隨后通過波峰,平臺前端下俯,縱搖由負轉(zhuǎn)正,整個過程縱搖曲線呈三角函數(shù)趨勢。由于尾部子平臺最后受到波浪沖擊,波浪能量衰減,因而縱搖振幅最小,縱搖角度變化幅度隨著子平臺數(shù)量增加而逐段遞減。

由表2可知:子平臺越靠近尾部,平臺的平均縱蕩幅值越大,平均縱搖幅值越小;由于第1段子平臺前端系泊于造波池底,運動范圍受到系泊纜繩的限制,因而平均垂蕩幅值的最大值出現(xiàn)于第2段子平臺。隨著子平臺數(shù)量的增加,縱蕩幅值逐漸降低,整體平臺更加穩(wěn)定。

表2 浮體平臺在一個周期內(nèi)的平均運動響應

4 結 論

本文利用無網(wǎng)格法計算波浪運動大變形的優(yōu)勢和離散元法模擬多體運動的特長,將無網(wǎng)格MPS法和離散元DEM法進行耦合,建立了流固(多體)耦合數(shù)值模擬方法;采用反演式造波形式,提升了造波水池的造波精度和穩(wěn)定性;利用阻尼式消波的技術,縮減了計算區(qū)域,提升了計算效率;基于系泊浮體數(shù)值模型,設置并實現(xiàn)了首尾系泊的多段分段浮體平臺在規(guī)則波沖擊下的運動的數(shù)值模擬。

研究結果表明,分段浮體平臺的子平臺位置越靠近整體平臺下游,垂蕩和縱搖響應逐漸削弱,縱蕩響應成為主導特征。大型分段浮體平臺在設計時可通過增加系泊纜繩等措施來降低前部子平臺的位移,進而減小前部子平臺的垂蕩和縱搖;其次,還應增加連接部件對尾部子平臺的約束力,防止極端海況下尾部子平臺間的相互碰撞以及尾部子平臺與主體平臺的脫離;同時,適當增加子平臺數(shù)量有利于減小整體平臺的縱蕩幅值,保持穩(wěn)定性。研究結果可為大型浮體平臺的設計及優(yōu)化提供參考和依據(jù)。