劉 圓， 張普然， 劉 璐， 季順迎*

(1.中國船級社海洋工程技術(shù)中心,天津 300457; 2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室,大連 116024)

1 引 言

在自然界與人類的生產(chǎn)活動中，顆粒-流體耦合現(xiàn)象廣泛存在[1]，如滑坡泥石流[2]、水下拋石及冰水耦合等，在地質(zhì)災(zāi)害演變、化工過程和巖土工程中具有重要的研究意義與價值。離散元方法(DEM)是模擬不同形態(tài)顆粒的主要方法，能夠清晰地描述顆粒離散介質(zhì)的動力學(xué)行為[3]。同時，求解流體的方法多樣，目前應(yīng)用較多的有有限元方法[4]、無網(wǎng)格粒子法[5]和格子-玻爾茲曼方法[6]等。采用數(shù)值模擬開展流固耦合計算是重要的研究手段[7]。

早期離散元法的單元形式主要是二維圓盤與三維圓球[8]，為了更好地描述形態(tài)各異的真實(shí)顆粒的力學(xué)行為，發(fā)展了許多非球形的幾何形態(tài)，如多面體[9]、超二次曲面[10]與擴(kuò)展多面體等。擴(kuò)展多面體單元能夠較好地模擬海冰和碎石等結(jié)構(gòu)。擴(kuò)展多面體結(jié)合了球體與多面體的幾何特征，避免了復(fù)雜且不穩(wěn)定的幾何判斷，提高了單元接觸搜索的效率和穩(wěn)定性[11]。由于單元表面為曲面，擴(kuò)展多面體間的接觸力計算可采用赫茲接觸模型[12]。此外，為分析工程中連續(xù)的脆性材料破碎為顆粒材料的過程，發(fā)展了基于擴(kuò)展多面體的粘結(jié)與破碎模型[13,14]?；跀U(kuò)展多面體的離散元方法可與多種計算流體的方法耦合，用于顆粒-流體耦合系統(tǒng)的模擬。

流固耦合中經(jīng)常存在劇烈沖擊，流體發(fā)生大變形。SPH方法采用無網(wǎng)格的拉格朗日粒子模擬流體的動力學(xué)過程，能夠有效避免網(wǎng)格畸變并可自適應(yīng)地模擬自由表面[15]，SPH方法和DEM方法的基本計算單元均是相似的離散粒子，數(shù)據(jù)結(jié)構(gòu)和算法邏輯具有較高的相似度，在程序編譯與計算方面具有顯著優(yōu)勢[16]。

DEM-SPH耦合模型已經(jīng)成功應(yīng)用于多個領(lǐng)域，包括海洋工程、生物力學(xué)與水利工程等。具體實(shí)例包括海浪沖擊消波塊[17,18]、植物細(xì)胞干燥過程[19,20]和毛細(xì)血管中血細(xì)胞的運(yùn)動[21]等。DEM-SPH方法的有效性和準(zhǔn)確性也在與試驗(yàn)和傳統(tǒng)數(shù)值結(jié)果的比較中得到了充分的驗(yàn)證[22-24]。

水下拋石過程是海洋工程和水利工程中常見的問題，拋石法是工業(yè)界保護(hù)水底管道的主要方法[25]。本文采用基于擴(kuò)展多面體的離散元方法模擬碎石塊，采用弱可壓縮的SPH方法模擬流體，并采用邊界排斥力模型計算擴(kuò)展多面體與SPH粒子間的相互作用，從而建立了碎石塊與水相互作用的DEM-SPH耦合方法。為模擬卸料斗中落石與水底管線的碰撞過程，建立了擴(kuò)展多面體和SPH粒子與錐形結(jié)構(gòu)的接觸搜索方法。采用DEM-SPH方法分析了碎石塊通過錐形結(jié)構(gòu)入水并與水底管道的作用過程，分析了卸料斗靜止與運(yùn)動時對管線的沖擊力和管線周圍落石堆積形態(tài)的影響。

2 DEM-SPH耦合模型

2.1 擴(kuò)展多面體單元模型

如圖1所示，對球體和任意形狀多面體進(jìn)行閔可夫斯基求和運(yùn)算，即可構(gòu)成擴(kuò)展多面體。其中，多面體稱為基本多面體，球體稱為擴(kuò)展球體，球體半徑稱為擴(kuò)展半徑。擴(kuò)展多面體將角點(diǎn)和棱邊的接觸搜索轉(zhuǎn)化為相應(yīng)的球面和圓柱面搜索，簡化了單元間的接觸判斷過程[26]。

圖1 基于Minkowski Sum方法的擴(kuò)展多面體單元

(1)

式中ri和rj為兩個擴(kuò)展多面體單元的擴(kuò)展半徑。

δi j<0時表示兩個擴(kuò)展多面體發(fā)生接觸，接觸力的計算采用非線性模型，包括法向力和切向力。法向力包含彈性力與粘滯力，可表示為

(2)

(3，4)

式中ks為單元間切向有效剛度，且ks=rs nkn，而rs n=1/(2+2ν)，ν為泊松比，δs為切向位移，μ為滑動摩擦系數(shù)。

2.2 弱可壓縮SPH方法

弱可壓縮SPH方法在流體動力學(xué)的數(shù)值模擬中得到廣泛應(yīng)用。流體粒子密度的連續(xù)性方程為

(5)

式中ui j為粒子i與j之間的速度差，ui j=ui-uj，W為核函數(shù)，本文采用三次樣條核函數(shù)。每個粒子的壓力可通過狀態(tài)方程計算，即

(6)

(7)

2.3 離散元顆粒與SPH粒子間的排斥力模型

SPH方法中常用的邊界粒子法和虛擬粒子法不適用于構(gòu)造復(fù)雜形狀的邊界[27]。因此，如圖2所示，用排斥力方法簡化了SPH粒子與擴(kuò)展多面體之間的相互作用，作用力只與距離d相關(guān)。

圖2 SPH粒子與固體邊界之間的排斥力模型

fB=βnBε(z)R(d)

(8)

式中fB為邊界排斥力，β為控制系數(shù)，nB為壁面邊界法向量，ε(z)為與粒子的z方向坐標(biāo)有關(guān)的壓力校正項(xiàng)，d為SPH粒子與邊界之間的最小距離，R(d)為排斥力函數(shù)，可表示為

(9)

2.4 擴(kuò)展多面體與錐形結(jié)構(gòu)間的接觸判斷

本文采用錐形結(jié)構(gòu)來模擬海底管線以及卸料斗，建立了錐形結(jié)構(gòu)與擴(kuò)展多面體的接觸判斷方法。由于擴(kuò)展多面體單元主要由角點(diǎn)、棱邊和平面等元素組成，錐形結(jié)構(gòu)外部與擴(kuò)展多面體單元的接觸可分為以下三種類型。(1) 球體與錐形結(jié)構(gòu)接觸；(2) 圓柱體與錐形結(jié)構(gòu)接觸；(3) 平面與錐形結(jié)構(gòu)接觸。錐形結(jié)構(gòu)內(nèi)部與擴(kuò)展多面體單元的接觸為球體與錐形結(jié)構(gòu)接觸。

圖3 球體與錐形結(jié)構(gòu)間的接觸

圖4 圓柱與錐形結(jié)構(gòu)的接觸判斷

圖5 平面與錐形結(jié)構(gòu)的接觸判斷

3 DEM-SPH耦合方法的驗(yàn)證

3.1 潰壩模擬

三維潰壩流的計算模型如圖6所示，水壩的z方向長度H=0.6 m，x方向長度L=2H，y方向長度W=3H。長方體容器的z方向長度D=3H，x方向長度d=5.366H，y方向長度與潰壩相同。SPH粒子光滑長度為0.015 m，粒子總數(shù)為25600。

圖6 潰壩計算模型

圖7 潰壩過程的SPH模擬

圖8 潰壩流前緣位置隨時間的變化

3.2 楔形塊入水模擬

楔形塊入水是結(jié)構(gòu)物入水問題的經(jīng)典算例，已有大量的工作對其開展了研究，包括理論、試驗(yàn)與數(shù)值模擬，因此通過該算例驗(yàn)證DEM-SPH耦合方法的合理性。參考Lewis等[30]的試驗(yàn)，模型尺寸如圖9所示，楔形塊的質(zhì)量為23.4 kg，平均密度為336.89 kg/m3。邊界長度為2.0 m，寬度比楔形塊寬度稍大一些，為0.76 m。水深0.59 m，楔形塊最底部位置距水面高度為0.5 m。

圖9 楔形塊入水模型

研究不同SPH粒子光滑長度下楔形塊加速度隨時間的變化，并與試驗(yàn)結(jié)果進(jìn)行了對比，如圖10所示。由于試驗(yàn)中存在空氣阻力與試驗(yàn)裝置的摩擦，楔形塊未入水之前的加速度略小于重力加速度。當(dāng)光滑長度取0.01 m或0.02 m時，模擬結(jié)果與試驗(yàn)結(jié)果較為接近，驗(yàn)證了DEM-SPH耦合方法的合理性。當(dāng)光滑長度擴(kuò)大到0.05 m時，加速度峰值相較于試驗(yàn)結(jié)果有明顯下降。實(shí)際上，DEM單元與SPH粒子的粒徑比越大，模擬結(jié)果越接近實(shí)際情況，但是SPH粒子過多時計算效率會顯著下降，因此DEM單元與SPH粒子的粒徑比一般保持在8～15之間。

圖10 不同光滑長度下的加速度

4 水底管道拋石加固的流固耦合模擬

鋪設(shè)在海床上的海底管道或電纜通常需碎石掩埋，以增強(qiáng)管道的穩(wěn)定性并保護(hù)管道，使其免受侵蝕與沖擊。該過程通常是通過船舶上的落石裝置將碎石傾倒在管道上[31]，如圖11所示。管道周圍落石的數(shù)量與堆積形態(tài)對于評估落石是否對管道提供了足夠的約束具有重要意義。采用本文提出的DEM-SPH耦合方法模擬落石過程，并研究卸料斗靜止與運(yùn)動時落石的堆積形態(tài)。

圖11 水底管線落石過程

4.1 模型建立與參數(shù)設(shè)置

采用Voronoi法(VTA)生成隨機(jī)形態(tài)的擴(kuò)展多面體用于模擬落石入水過程。根據(jù)3.2節(jié)的結(jié)論，流固耦合模擬過程中DEM單元與SPH粒子的粒徑比一般保持在8～15之間。綜合考慮實(shí)際的結(jié)構(gòu)參數(shù)與計算效率，確定模擬過程中SPH粒子的光滑長度為0.05 m，DEM單元的粒徑約為 0.4 m～0.6 m。

如圖12所示，卸料斗由自上而下的三段錐形結(jié)構(gòu)組成。最下面的落石管深入水面下方，頂部管道用于容納落石在卸料斗中的堆積。水底管道的直徑為0.8 m，其軸線方向?yàn)閥方向。卸料斗和水底管道都是剛性結(jié)構(gòu)。所有石塊最初都放置在卸料斗的頂部管道中，初始狀態(tài)有豎直向下的初速度，并在重力作用下從管道中掉落。

圖12 卸料斗

表1 落石模擬中的主要參數(shù)

4.2 卸料斗固定的模擬結(jié)果與分析

卸料斗的三根錐形結(jié)構(gòu)不發(fā)生移動，底部水域邊界為4 m×4 m，水深為3 m。該模擬用于驗(yàn)證DEM單元以及SPH粒子與錐形結(jié)構(gòu)間的接觸模型。圖13顯示了落石的模擬過程。其中，石塊下落初速度為1 m/s，與卸料斗內(nèi)部發(fā)生接觸并落入水中，然后與水底管道發(fā)生碰撞，最終落石完全覆蓋水底管道。隨著水中石塊數(shù)量的增多，卸料斗底部產(chǎn)生堵塞現(xiàn)象，石塊在卸料斗中產(chǎn)生了堆積現(xiàn)象。

圖13 管線落石的流固耦合模擬過程

圖14顯示了水底管道受到的沖擊力，包括落石和水對管道的沖擊力。開始一段時間管道只受到水的沖擊力，1.2 s開始有落石與管道發(fā)生接觸，由于邊界較為狹小，在t=13 s之后，落石在卸料斗中堆積，無法入水，水底管道附近的落石沒有繼續(xù)增多，因此落石對管的沖擊力并未隨時間一直增長，反而在堆積穩(wěn)定后出現(xiàn)下降趨勢。開始0～2 s內(nèi)水對管道的沖擊力有一個較大的波動，穩(wěn)定后水對管道的力在25 kN左右。整個過程中，落石對管道的力基本處于10 kN～30 kN。

圖14 水底管道受到的沖擊力

4.3 卸料斗移動的模擬結(jié)果與分析

通過4.2節(jié)卸料斗固定的模擬，成功驗(yàn)證了石塊以及水與卸料斗間的接觸模型。而在實(shí)際情況中，拋石船會沿著水底管線不斷移動，卸料斗也隨船一起移動。結(jié)合實(shí)際情況與計算效率，確定水域大小等參數(shù)，模擬了卸料斗邊移動邊拋石的過程。

卸料斗沿著管道軸線方向移動。底部水域邊界尺寸為4 m×10 m，水深為3.0 m。在模擬過程中，落石下落初速度與卸料斗移動速度是兩個重要參數(shù)，影響落石最終的堆積形態(tài)。如圖15所示，當(dāng)卸料斗移動速度較小而落石初速度較大時，很容易出現(xiàn)落石堵塞卸料斗的情況。因此，選取合適的落石初速度與卸料斗移動速度，是成功模擬水下拋石過程的關(guān)鍵。

圖15 落石堵塞卸料斗

圖16顯示了卸料斗速度為0.5 m/s，落石初速度為1 m/s時的拋石過程。卸料斗開始停止，2.36 s后向y軸正方向做勻速直線運(yùn)動。石塊通過卸料斗落入水中，然后撞擊水底管道。水底管道上的作用力如圖17所示。隨著時間的增長，水對管道的作用力趨于穩(wěn)定，而落石對管道的作用力總體呈不斷上升趨勢，這表明石塊傾倒過程中未出現(xiàn)堵塞的情況。

圖16 卸料斗移速0.5 m/s時的模擬結(jié)果

圖17 卸料斗移速0.5 m/s時管道受力

如圖16(b)所示，拋石結(jié)束后水底管道并未完全受落石覆蓋，這種情況下落石無法充分保護(hù)管道，是工程中必須要避免的現(xiàn)象。保持卸料斗移動速度不變，增大落石初速度，觀察拋石結(jié)束后落石的分布情況。

圖18展示了落石下落初速度為2 m/s和 3 m/s 時的模擬結(jié)果?？梢钥闯?，在卸料斗移速為0.5 m/s的條件下，落石初速度大于2 m/s時，落石可以完全覆蓋管道，15 s后卸料斗到達(dá)邊界位置，不再移動，此時落石已經(jīng)開始在卸料斗中堆積，落石對水底管道的沖擊力不再增大，與4.2節(jié)的情況相同。

圖18 改變下落初速度后的模擬結(jié)果

當(dāng)落石初速度取3 m/s和2 m/s時，拋石結(jié)束后管道所受沖擊力比較接近。這表明在落石已經(jīng)能夠覆蓋管道時，隨著落石初速度增大，沖擊力增大的幅度在減小。由于受拋石模型中邊界約束以及卸料斗底部到水底管道的距離限制，覆蓋管道的落石數(shù)量存在上限。而拋石結(jié)束后，管道所受沖擊力主要是穩(wěn)定后覆蓋管道落石的重力，因此可以預(yù)測，對于本文給出的拋石模型，落石初速度大于 3 m/s 后，落石在管道附近的分布形態(tài)以及管道所受沖擊力的變化會越來越小。

圖19 落石初速度不同時落石對管道的沖擊力

接下來討論卸料斗移動速度對模擬過程的影響。落石下落初速度固定為3 m/s，圖20展示了卸料斗移動速度為1 m/s和1.5 m/s時的模擬結(jié)果。相較于圖18(b)，圖20中覆蓋管道的落石減少很多，圖20(b)所示的結(jié)果已經(jīng)屬于拋石失敗的范疇。

圖20 卸料斗移動速度不同時的模擬結(jié)果

由于落石初速度相同，單位時間內(nèi)下落的落石數(shù)量相同，因此拋石過程中管道受力隨時間的變化曲線基本重合，因此將自變量更改為卸料斗移動的距離。分析了落石對管道沖擊力隨卸料斗移動距離的變化。如圖21所示，卸料斗開始移動時，三個工況下沖擊力很接近，拋石過程中，卸料斗移動速度越快，拋石的總量越小，管道受沖擊力也越小。其中，卸料斗移速為1 m/s和1.5 m/s時沖擊力比較接近，是由于這兩種工況下，覆蓋在管道上方的落石數(shù)量太少，而散落在管道兩側(cè)的落石對管道的沖擊力基本可以忽略不計。

圖21 卸料斗移動速度不同時落石對管道的沖擊力

對于本文構(gòu)建的水下拋石模型，通過幾組工況的模擬，可以確定當(dāng)落石初速度取2 m/s～3 m/s，卸料斗移動速度取0.5 m/s～1 m/s時，卸料斗移動拋石結(jié)束后落石分布狀況較為理想，可以實(shí)現(xiàn)管道保護(hù)的目標(biāo)。落石初速度過大會造成堵塞，過小會導(dǎo)致落石無法保護(hù)管道，而卸料斗移速對模擬結(jié)果的影響與落石初速度相反。在實(shí)際的拋石過程中，可以先通過本文提出的卸料斗拋石模型進(jìn)行多工況的數(shù)值模擬，根據(jù)結(jié)果確定相關(guān)參數(shù)后再進(jìn)行作業(yè)，可以有效節(jié)省人力物力。

5 結(jié) 論

本文采用基于閔科夫斯基的擴(kuò)展多面體單元建立隨機(jī)形態(tài)顆粒的離散元方法，采用弱可壓縮格式的SPH方法進(jìn)行水動力學(xué)模擬，并采用排斥力模型計算擴(kuò)展多面體與SPH粒子之間的相互作用力，從而建立了DEM-SPH耦合模型。通過潰壩及楔形塊入水等典型算例驗(yàn)證了本文方法的穩(wěn)定性和合理性。采用DEM-SPH耦合方法模擬了水底管道落石過程，對水底管道受到的作用力和碎石塊堆積形態(tài)進(jìn)行了分析。設(shè)定多組工況，探究了落石初速度和拋石船移動速度對拋石結(jié)果的影響。計算結(jié)果表明，本文建立的DEM-SPH方法可有效模擬落石入水并與水底管道作用的過程，可進(jìn)一步應(yīng)用于落石裝置設(shè)計和水底管道穩(wěn)定性等工程問題的分析中。