基于OpenFOAM與PFC耦合方法的水下滑坡數(shù)值模擬研究*

陳聞瀟 石 崇 單治鋼 張一平 李汪洋

(①河海大學(xué)， 巖石力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室， 南京 210098， 中國(guó)) (②河海大學(xué)， 巖土工程科學(xué)研究所， 南京 210098， 中國(guó)) (③中國(guó)電建集團(tuán)， 華東勘測(cè)設(shè)計(jì)研究院有限公司， 杭州 310014， 中國(guó))

0 引 言

水下滑坡作為在近陸范圍常見(jiàn)的一種地質(zhì)災(zāi)害，近年來(lái)已造成嚴(yán)重的人員傷亡或財(cái)產(chǎn)損失(Alves et al.，2015; Terry et al.，2017;Sassa et al.，2019)，僅在1616～1886年期間，世界上由于地震、火山噴發(fā)等引發(fā)的水下滑坡失穩(wěn)事件達(dá)333次(Milne，1897)。此外水下滑坡對(duì)海洋設(shè)施如海洋油氣生產(chǎn)平臺(tái)、水下管道、水下電纜等也具有較大威脅(朱超祁等， 2016)，因此研究水下滑坡的形成過(guò)程與機(jī)理具有重要的工程意義。

由于水下滑坡的特殊環(huán)境，數(shù)值模擬是研究水下滑坡的重要手段。楊林青(2012)基于ABAQUS有限元強(qiáng)度折減法分析了水下斜坡穩(wěn)定性，霍沿東等(2019)基于極限分析上限方法分析了水下斜坡在地震下的穩(wěn)定性。但這些方法采用的連續(xù)力學(xué)分析方法，無(wú)法模擬分析滑坡造成的大變形破壞、滑坡運(yùn)動(dòng)過(guò)程、堆積形態(tài)與災(zāi)害影響范圍等(陳聞瀟等， 2020)。

顆粒離散元(PFC)方法具有模擬大變形破壞的優(yōu)勢(shì)而被廣泛運(yùn)用到滑坡過(guò)程分析中(陳曉等， 2020)，但目前基于PFC方法的滑坡研究大多為陸上滑坡(曹文等， 2017； 周禮等， 2019)，對(duì)水下滑坡研究不足。且在PFC中無(wú)法直接考慮水的作用，常見(jiàn)的達(dá)西滲流法(畢慶杰， 2018)、管道域法(Zhang et al.，2019;Zhang et al.，2020)分別適用于多孔介質(zhì)滲流和水力劈裂，但無(wú)法模擬水下大尺度流體環(huán)境，故需要考慮引入流體力學(xué)計(jì)算方法，才能實(shí)現(xiàn)流固耦合計(jì)算。

OpenFOAM是一款計(jì)算流體力學(xué)的開(kāi)源軟件，相比商業(yè)軟件，其具有良好的擴(kuò)展性、豐富的模型庫(kù)和求解高效的優(yōu)勢(shì)(Shojaeefard et al.，2017)。已經(jīng)用來(lái)研究波浪、風(fēng)場(chǎng)等方面的流體作用(Jacobsen et al.，2015; Liu et al.，2017)，但其在土木和地質(zhì)方面鮮有研究，有待進(jìn)一步開(kāi)發(fā)。

本文通過(guò)PFC內(nèi)置TCP套接字模塊(Itasca， 2014)，實(shí)現(xiàn)基于OpenFOAM與PFC耦合計(jì)算方法，并針對(duì)PFC流體分析提出適用于該耦合方法的相似性方法。然后通過(guò)典型案例，驗(yàn)證該耦合方法進(jìn)行水下滑坡分析的可行性，進(jìn)一步研究水下滑坡的動(dòng)力學(xué)特性及堆積形態(tài)，并與單向耦合水下滑坡和陸上滑坡結(jié)果進(jìn)行對(duì)比分析。

1 耦合分析實(shí)現(xiàn)方法

在考慮滑坡與水動(dòng)力學(xué)耦合作用時(shí)，流體動(dòng)力學(xué)需要考慮連續(xù)性方程和動(dòng)量守恒方程，此處采用考慮顆粒作用修正后的連續(xù)性方程和Navier-Stokes方程(Jackson， 2000)：

(1)

(2)

式中：V*為流體速度矢量；n為流體單元內(nèi)的孔隙率;ρf為流體密度；p為水壓力；t為時(shí)間；τ為黏性力張量；g為重力加速度；fint為單個(gè)流體網(wǎng)格內(nèi)包含顆粒所受的拖曳合力。上述方程的計(jì)算使用OpenFOAM中的demIcoFoam求解器，demIcoFoam通過(guò)修改icoFoam求解器來(lái)模擬固體顆粒的作用，同時(shí)添加顆粒與流體之間的相互作用力和基于粗網(wǎng)格法(Wang et al.，2013)得到的單元孔隙率。

單個(gè)體積為Vi的網(wǎng)格內(nèi)的拖拽合力fint通過(guò)PFC中的CFD模塊計(jì)算，計(jì)算公式為：

(3)

式中：fdrag為單個(gè)顆粒所受流體對(duì)其的拖曳力，其通過(guò)下式計(jì)算：

(4)

式中：r為顆粒半徑；v為流體速度；u為顆粒速度?？紫堵蕁通過(guò)基于粗網(wǎng)格法的多面體法計(jì)算，n-χ項(xiàng)是考慮局部孔隙度的經(jīng)驗(yàn)系數(shù)。這個(gè)修正項(xiàng)使拖曳力同時(shí)適用于高孔隙和低孔隙度系統(tǒng)，流體雷諾數(shù)也可在大范圍內(nèi)取值。拖拽力系數(shù)Cd被定義為：

(5)

式中：當(dāng)流體的動(dòng)力黏滯系數(shù)為μf時(shí)，顆粒的雷諾數(shù)Rep通過(guò)下式計(jì)算：

(6)

上述耦合方法，基于PFC中提供的TCP套接字(Itasca， 2014)連接OpenFOAM與PFC，實(shí)現(xiàn)耦合計(jì)算中信息傳遞與更新。每一個(gè)流程步中，OpenFOAM將流體速度v、壓強(qiáng)p和壓力梯度p傳遞給PFC，而后在PFC中進(jìn)行式(4)～式(6)的計(jì)算以得出每個(gè)顆粒所受拖拽力fdrag，繼而通過(guò)力-位移運(yùn)動(dòng)方程得出土體顆粒的運(yùn)動(dòng)狀態(tài)，并通過(guò)式(3)～式(6)更新單元格內(nèi)拖拽合力fint，將其與單元格孔隙率n一起傳遞給OpenFOAM實(shí)現(xiàn)式(1)～式(2)的計(jì)算，這樣便完成了一個(gè)耦合循環(huán)步的計(jì)算。耦合流程如圖 1 所示：

圖 1 耦合交互方法流程圖Fig. 1 Flowchart of coupling interaction method

用顆粒水下自由沉降來(lái)驗(yàn)證該方法的可行性(Jing et al.，2019)，將半徑為5 mm的土體顆粒置于流體環(huán)境中，計(jì)算顆粒在重力下自由運(yùn)動(dòng)0.5 s后的運(yùn)動(dòng)速度。圖 2為運(yùn)動(dòng)速度隨時(shí)間的變化曲線(xiàn)，其理論解由式(7)計(jì)算，可得顆粒速度的理論解與模擬值一致，可驗(yàn)證該流固耦合方法正確性。

(7)

圖 2 不同時(shí)刻的顆粒下落速度Fig. 2 Falling velocities of particles at different times

在流固耦合中由于模型尺寸或顆粒尺寸相比真實(shí)值變化，雷諾數(shù)、流速、壓力差、啟動(dòng)力等均會(huì)變化，需要考慮相似性理論(孫其誠(chéng)等， 2008； 蔣明鏡等， 2014)?，F(xiàn)有的相似性方法主要有兩種：粒徑縮放-流體參數(shù)調(diào)整法和粒徑恒定-流體參數(shù)調(diào)整法(倪小東等， 2009)。在離散元計(jì)算中，通常以大的離散顆粒來(lái)替代真實(shí)試驗(yàn)中的土體顆粒(Liu et al.，2015)，故相當(dāng)于增大了顆粒尺寸同時(shí)模型尺寸不變，因此現(xiàn)有的相似性方法無(wú)法很好地模擬這種情況。通過(guò)調(diào)節(jié)流體參數(shù)的方法均無(wú)法較好地模擬，其會(huì)造成流場(chǎng)結(jié)果失真，故需要提出一種基于離散元PFC的相似性方法。

當(dāng)顆粒尺寸為真實(shí)值的N倍時(shí)，本文采取程序二次開(kāi)發(fā)的方式，通過(guò)修改內(nèi)部函數(shù)使得模型滿(mǎn)足相似性原理。首先通過(guò)二次開(kāi)發(fā)，將雷諾數(shù)Rep賦以系數(shù)1/N，用以保證顆粒半徑擴(kuò)大N倍的同時(shí)雷諾數(shù)(式6)不變，即滿(mǎn)足雷諾準(zhǔn)則。同時(shí)將拖曳力fdrag賦以系數(shù)N，用以保證拖拽力與重力同步增長(zhǎng)，滿(mǎn)足啟動(dòng)力相似(倪小東等， 2009)，同時(shí)整體滿(mǎn)足弗洛德(Froude)準(zhǔn)則v2/gl=c1和歐拉準(zhǔn)則Δp/ρv2=c2。其中拖曳力fdrag的擴(kuò)大通過(guò)PFC中FISH函數(shù)二次開(kāi)發(fā)，在CFD_AFTER_UPDATE注冊(cè)點(diǎn)修改函數(shù)值實(shí)現(xiàn)。此外保持其他流體參數(shù)不變。

表 1顯示了不同半徑的顆粒，在通過(guò)該相似性方法模擬后，顆粒在重力下自由運(yùn)動(dòng)0.5 s后的平均速度和平均位移。為保證放大倍數(shù)變量的單一性，防止孔隙率變化對(duì)結(jié)果的干擾，故取顆?？傮w積一致。由表 1 可見(jiàn)顆粒平均速度和平均位移基本一致，誤差較小，可認(rèn)為該相似性方法具有一定的可行性與準(zhǔn)確性，并且通過(guò)放大顆粒半徑，可以縮小模擬所需顆粒個(gè)數(shù)，極大地縮短計(jì)算所需時(shí)間。

表 1 相似性參數(shù)及驗(yàn)證結(jié)果Table 1 Similarity parameters and results

2 水下滑坡的耦合計(jì)算分析

基于典型水下滑坡案例，分析OpenFOAM與PFC耦合水下滑坡模擬方法的可行性，并進(jìn)行水下滑坡的動(dòng)力學(xué)特性及堆積形態(tài)研究。

為兼顧問(wèn)題的三維條件和耦合方法的可拓展性，耦合方法建立于三維條件下，可以通過(guò)假三維模型來(lái)進(jìn)行平面問(wèn)題的模擬，如圖 3 所示建立連續(xù)邊坡模型，假三維厚度為1 m?；陔S機(jī)種子法生成水下滑坡的離散元模型，邊坡模型長(zhǎng)120 m，高25 m，共生成顆粒6747個(gè)，顆粒粒徑大小在0.25～0.5 m之間，呈線(xiàn)性分布。邊坡右側(cè)長(zhǎng)80 m，使得滑坡體能夠充分地運(yùn)動(dòng)堆積，同時(shí)避免在流體計(jì)算時(shí)邊界效應(yīng)的影響。對(duì)水下土體的模擬采用接觸黏結(jié)模型，在顆粒破碎之后，接觸模型由接觸黏結(jié)模型，退化成線(xiàn)性接觸模型，能夠較好地模擬土體的力學(xué)響應(yīng)和破壞特性(Shimizu， 2004)。本文模型參數(shù)如表 2 所示。

圖 3 初始計(jì)算模型及網(wǎng)格劃分Fig. 3 Initial calculation model and mesh generation

表 2 計(jì)算模型細(xì)觀力學(xué)參數(shù)Table 2 The meso mechanical parameters of calculation model

采用固定粗糙網(wǎng)格法建立流體網(wǎng)格模型，網(wǎng)格尺寸大小為1 m 的正方形，單元數(shù)量為7200。同時(shí)假定研究區(qū)域內(nèi)水下土體平均粒徑為0.001～0.002 m之間，結(jié)合上文提出的放大系數(shù)為250，在OpenFOAM的demIcoFoam求解器中設(shè)置流體的密度ρf=1000 kg · m-3，動(dòng)力黏滯系數(shù)μf=0.001 Pa · s。此外在耦合模塊中，在CFD_AFTER_UPDATE注冊(cè)點(diǎn)對(duì)雷諾數(shù)Rep和拖曳力fdrag賦以放大系數(shù)。

水下滑坡的耦合計(jì)算結(jié)果如圖 4 所示，后文描述時(shí)間為滑坡實(shí)際時(shí)間。圖 4a，圖4c分別為模型計(jì)算至29.9 s和39.8 s時(shí)的土體位移圖，可見(jiàn)滑坡呈現(xiàn)出較為明顯的圓弧滑動(dòng)面，隨著滑坡的發(fā)展，滑坡土體沿滑動(dòng)面運(yùn)動(dòng)，滑動(dòng)區(qū)域下部的土體不斷推進(jìn)，邊坡后緣土體逐漸破壞，土體最大位移分別為26.0 m和33.9 m，其堆積至坡腳的范圍分別為14.3 m和20.1 m。值得注意的是，滑坡體前端厚度較大，且前端呈現(xiàn)橢圓狀，同時(shí)在滑坡體的前端底部出現(xiàn)了滑水現(xiàn)象(修宗祥等， 2016)，其主要由于滑坡體快速運(yùn)動(dòng)而在底部形成了類(lèi)似水楔的水層，這同時(shí)也和景路等(2019)研究中得到的現(xiàn)象相符。由圖 4b 可見(jiàn)滑坡土體運(yùn)動(dòng)帶動(dòng)流體運(yùn)動(dòng)并形成了漩渦狀，流場(chǎng)速度分布基本與土體所在區(qū)域重合，遠(yuǎn)離滑坡區(qū)域的流體速度較小。隨著滑坡的發(fā)展，由圖 4d可見(jiàn)邊坡后緣部分由于土體滑落，原本由土體填充的區(qū)域由流體涌入，并在該區(qū)域內(nèi)形成了雜亂的紊流。

圖 4 耦合模型計(jì)算結(jié)果Fig. 4 The calculation results of the coupling modela. 29.9 s位移圖(單位：m)； b. 29.9 s流體速度矢量圖(單位：m · s-1)； c. 39.8 s位移圖(單位：m)； d. 39.8 s流體速度矢量圖(單位：m · s-1)

滑坡結(jié)束后堆積狀態(tài)如圖 5a 所示，滑坡體最終在坡腳堆積并停止進(jìn)一步運(yùn)動(dòng)，堆積距離為27.6 m，堆積體前端變薄，不再呈現(xiàn)橢圓面。由圖 5b 可見(jiàn)滑坡體孔隙率要明顯高于未滑動(dòng)土體，其孔隙率值處于0.4～0.5之間。這說(shuō)明隨著滑坡的發(fā)展，土體由緊密接觸的整體，破碎成了松散的團(tuán)體。

圖 5 滑坡結(jié)束后堆積狀態(tài)Fig. 5 The accumulation state after the landslidea. 滑坡體位移圖(單位：m)； b. 網(wǎng)格孔隙率圖

滑動(dòng)土體的速度位移曲線(xiàn)如圖 6 所示，可見(jiàn)在滑坡發(fā)展初期速度迅速增大，在15 s左右土體的平均速度達(dá)到峰值0.47 m · s-1，峰值的迅速出現(xiàn)是因?yàn)榱黧w對(duì)顆粒的拖拽力fdrag與顆粒速度方向相反，大小與流固速度差相關(guān)，故會(huì)限制土體速度進(jìn)一步發(fā)展，同時(shí)由于流固之間的相互作用，滑坡土體的能量一部分消耗于帶動(dòng)流體的運(yùn)動(dòng)。在20～80 s內(nèi)土體速度不斷減小，在50 s時(shí)土體平均速度下降到0.091 m · s-1左右，之后土體進(jìn)入緩慢的堆積平衡階段，此時(shí)流體對(duì)顆粒的拖拽力fdrag會(huì)減緩?fù)馏w顆粒速度下降的趨勢(shì)。如果該研究對(duì)象處于有海洋洋流的存在情況下，將會(huì)使得堆積體在流體的作用下運(yùn)移出較遠(yuǎn)的距離，這將會(huì)導(dǎo)致滑坡影響范圍更進(jìn)一步地?cái)U(kuò)大。

圖 6 滑坡體的平均速度位移曲線(xiàn)Fig. 6 Average speed and displacement curve of landslide body

3 耦合對(duì)比分析

為對(duì)水下滑坡的耦合計(jì)算方法進(jìn)一步探究，并討論水下滑坡與陸上滑坡的異同點(diǎn)，本節(jié)研究基于OpenFOAM-PFC的雙向耦合水下滑坡、單向耦合水下滑坡與陸上滑坡(無(wú)流固耦合)3種不同的耦合方式。單向耦合是當(dāng)前滑坡分析常用方法之一(丁歡歡等， 2013)，其在計(jì)算過(guò)程中信息只由流場(chǎng)傳遞至顆粒速度場(chǎng)，而顆粒運(yùn)動(dòng)不會(huì)對(duì)流體產(chǎn)生影響，即后續(xù)計(jì)算中流體速度v恒定不變(式4)。兩個(gè)模型均采用上文的模型尺寸與土體參數(shù)，但由于陸上滑坡初始條件下不受水力作用，故需要通過(guò)強(qiáng)度折減法調(diào)整邊坡的接觸強(qiáng)度而導(dǎo)致其產(chǎn)生滑坡。

可見(jiàn)陸上滑坡( 圖 7a) 和流固單向耦合滑坡( 圖 7b) 與基于雙向耦合的水下滑坡( 圖 4) 滑坡?tīng)顟B(tài)并不一致，兩者在滑坡過(guò)程中滑坡體前端厚度較小，均沒(méi)有橢圓面出現(xiàn)，無(wú)滑水現(xiàn)象產(chǎn)生，這是因?yàn)閱蜗蝰詈夏M中顆粒對(duì)流場(chǎng)作用無(wú)法體現(xiàn)，流固相互作用機(jī)理無(wú)法實(shí)現(xiàn)。陸上滑坡由于沒(méi)有流體作用，滑坡時(shí)間更短，主要位移發(fā)生在0～10 s內(nèi)，滑坡過(guò)程中的顆粒碰撞相對(duì)劇烈，大部分接觸黏結(jié)已經(jīng)破壞，滑坡體相對(duì)更加松散。而單向耦合中由于流體速度恒定不變，較大的流固流速差使得土體運(yùn)動(dòng)緩慢，滑坡過(guò)程持續(xù)達(dá)300 s，緩慢的運(yùn)動(dòng)過(guò)程使碰撞和剪切相對(duì)較弱，土體相對(duì)仍較密實(shí)。

圖 7 滑坡位移曲線(xiàn)與運(yùn)動(dòng)狀態(tài)圖Fig. 7 Diagram of landslide displacement curve and motion statea. 陸上滑坡； b. 流固單向耦合

從圖 7a 可以看出，不同耦合方式模擬的運(yùn)動(dòng)規(guī)律并不一致。雙向耦合和陸上滑坡的平均速度峰值分別為0.47 m · s-1和1.95 m · s-1，出現(xiàn)的時(shí)間分別在14.3 s和3.0 s。而單向耦合中并不存在明顯的速度峰值階段，其滑坡速度在10～200 s內(nèi)維持在0.15～0.2 m · s-1范圍內(nèi)，并隨著時(shí)間逐漸降低。這是由于其假定流體速度v恒定不變，隨著顆粒速度u增大，拖拽力fdrag不斷增大進(jìn)而當(dāng)合外力趨于0時(shí)，顆粒速度不再明顯增加，且只因顆粒接觸力變化而小幅度變化，以上說(shuō)明了單向耦合方法具有一定的局限性。受流體作用影響，水下滑坡的滑坡發(fā)展與運(yùn)移要明顯滯后于陸上滑坡。從圖 8b 可以看出，單向耦合模擬的位移要小于雙向耦合，這是由于單向耦合中更大的流固速度差，導(dǎo)致更多的能量被消耗在了流固相互作用中。而陸上滑坡由于運(yùn)動(dòng)劇烈，接觸破壞更多，顆粒碰撞更加劇烈，導(dǎo)致能量損耗迅速，故平均位移要小于雙向耦合結(jié)果?？傮w而言，陸上滑坡的位移要小于水下滑坡，且真實(shí)情況下由于水下存在水流作用，其可能會(huì)帶動(dòng)滑坡土體進(jìn)一步運(yùn)動(dòng)，產(chǎn)生更大的破壞范圍。

圖 8 滑坡體速度和位移曲線(xiàn)Fig. 8 Velocity and displacement curves of the landslidea. 平均速度曲線(xiàn)； b. 平均位移曲線(xiàn)

4 結(jié) 論

本文基于PFC和OpenFOAM建立了流固耦合計(jì)算方法，提出了適用于該耦合方法的相似性賦值方法，利用案例分析了水下滑坡的動(dòng)力學(xué)特性及堆積形態(tài)，并與基于單向耦合方法的水下滑坡和陸上滑坡結(jié)果進(jìn)行了對(duì)比，得到主要結(jié)論如下：

(1)顆粒水下自由沉降模擬得到的計(jì)算值與理論計(jì)算值一致，表明了該耦合方法的正確性，同時(shí)對(duì)不同粒徑的顆粒采用相似性方法考慮，驗(yàn)證了該方法的可行性。

(2)水下滑坡的運(yùn)動(dòng)形態(tài)表現(xiàn)為滑坡體前端厚度較大，呈橢圓面并有滑水現(xiàn)象的產(chǎn)生，由于流固相互作用，拖拽力會(huì)阻礙顆粒運(yùn)動(dòng)，滑坡土體運(yùn)動(dòng)速度較快達(dá)到峰值并緩慢下降，堆積土體孔隙率明顯大于未滑動(dòng)土體。結(jié)論表明該耦合方法能夠較好地進(jìn)行水下滑坡過(guò)程模擬。

(3)水下滑坡與陸上滑坡在運(yùn)動(dòng)特征、堆積形態(tài)和滑坡歷程上具有較大的差異。單向耦合方法存在一定局限性，PFC-OpenFOAM雙向耦合能更好模擬水下滑坡運(yùn)動(dòng)特征。