穩(wěn)定流作用下的二維海底管道底部沖刷數(shù)值分析

張龍珠，陳善群，廖 斌

(安徽工程大學(xué) 力學(xué)重點實驗室，安徽 蕪湖 241000)

隨著科學(xué)技術(shù)的推陳出新以及地表資源的不斷消耗，人們將目光逐步投向了覆蓋地球面積71%且蘊含豐富能源的海洋。海底管道在人們開采海洋資源的過程中扮演著重要角色，承擔起輸送石油、天然氣等能源的艱巨任務(wù)。隨著世界工業(yè)規(guī)模的不斷發(fā)展以及能源消耗的逐步提高，海底管道的鋪設(shè)數(shù)量越來越多。但同時，海底管道也面臨著海洋自然環(huán)境的潛在威脅。具體表現(xiàn)在:海洋環(huán)境中水流運動的復(fù)雜多變驅(qū)動海底管道周圍泥沙的運動，使得沙床表面的泥沙被卷起并持續(xù)輸運，在管道底部形成明顯的沖刷坑。隨著沖刷坑的不斷發(fā)展加深，海底管道逐漸下垂會引起應(yīng)力集中，致使管道的懸挑部分會因應(yīng)力過度或疲勞而發(fā)生損壞。顯然，研究水流作用下海底管道的底部沖刷具有重要的現(xiàn)實意義和學(xué)術(shù)價值。

目前已有較多關(guān)于水流作用下海底管道底部沖刷的研究成果見諸發(fā)表。實驗研究方面，學(xué)者們分別從沖刷開始的臨界條件、尾流沖刷、Keulegan-Carpenter(KC)數(shù)、Shields數(shù)、多管道以及海底與管道間隙等方面[1-6]對管道沖刷過程進行了實驗研究。這些研究工作雖然給出了一些管道沖刷坑的形成規(guī)律、影響因素及經(jīng)驗公式，但受到空間尺度的限制，與工程實際存在較大差距。數(shù)值研究方面，Nadeem Ahmad[7]等采用開源流體動力學(xué)計算軟件REEF3D，建立了同向波流作用下管道沖刷的流體動力學(xué)計算模型，研究了海底管道下方最大沖刷深度的變化規(guī)律，但并未更深入地分析其他參數(shù)對沖刷坑深度的影響。史舒婧[8-9]運用FLUENT自定義函數(shù)(UDF)模塊建立了研究海底管道沖刷的數(shù)值模型，著重研究了Shields數(shù)對沖刷坑發(fā)展以及流場變化的影響，但是由于其對空間、時間離散以及邊界條件的一些簡化設(shè)定，使得數(shù)值模擬精度受到一定影響。胡鑫煒[10]利用Open FOAM結(jié)合修正的k-ε湍流模型建立了兩相泥沙沖刷模型，模擬穩(wěn)定來流下海底管道的沖刷過程，主要探索了埋置深度對管道沖刷坑形態(tài)及深度的影響。此外，周麗丹[11]等采用數(shù)值模擬結(jié)合實驗研究的方法對均勻來流海底管道的沖刷機理以及管道所受振動的特性進行了深入研究，解釋了沖刷坑的形成過程，找出了最大沖刷深度所在位置規(guī)律，但其總體研究成果主要依據(jù)實驗研究，對數(shù)值方面的探討不夠。

研究擬建立一個求解精度較高，適用于海底管道沖刷的數(shù)值計算模型，在對該數(shù)值模型進行驗證的基礎(chǔ)上，分析管道底部沖刷坑初步形成及發(fā)育機理。為深入發(fā)掘管道底部沖刷的影響因素，對不同流速、水深、管道直徑及管道與沙床間隙等條件下，管道底部沖刷坑的位置、深度等規(guī)律進行系統(tǒng)研究。研究工作有望為揭示海底管道底部沖刷機理提供信息參考。

1 數(shù)值模型

1.1 控制方程

在笛卡爾坐標系下，通過求解不可壓縮流體的雷諾時均納維-斯托克斯(RANS)方程,結(jié)合RNGk-ε湍流模型計算獲得數(shù)值水槽中的流場?？刂品匠痰亩S形式如下:

連續(xù)性方程：

(1)

動量方程：

(2)

式中,Vf是水體的體積分數(shù)；ρ是水體密度，取值1 000 kg/cm3；p為壓強；t為時間；u,v是流體沿x,y方向的速度分量；Ax,Ay是x,y方向水體的面積分數(shù)；RSOR是質(zhì)量源項；R,ξ是修正系數(shù)，與坐標系的選擇有關(guān)；Gx,Gy是流體在x,y方向上的重力加速度；fx,fy是流體在x,y方向上的粘滯力加速度。

(3)

(4)

式中，P、G均為湍流動能項(P由水體流速梯度引起，G由浮力引起)；kT、εT為擴散項。

1.2 FAVOR處理復(fù)雜模型

為精準刻畫海底管道的曲面外形、沙床隨水流沖刷而產(chǎn)生的形態(tài)變化過程以及節(jié)約計算資源，選取FAVOR(Fractional Area/Volume Obstacle Representation)網(wǎng)格處理技術(shù)進行建模。FAVOR方法即通過單元面積/體積分數(shù)精確描述物體復(fù)雜外型的方法，是一種獨特的網(wǎng)格處理技術(shù)。FAVOR方法與有限差分法處理曲面模型的對比如圖1所示。從圖1中可見，處理相同的幾何模型，F(xiàn)AVOR往往僅需較少的網(wǎng)格就可以達到很高的精度要求，但是傳統(tǒng)的有限差分法處理曲面模型卻要求更多的網(wǎng)格數(shù)量，需要更多的計算時間。即使是曲面造型的模型，也可以用矩形網(wǎng)格準確刻畫幾何外型，確保分析模型不會失真。

圖1 FAVOR方法與有限差分法處理曲面模型的對比

FAVOR方法在求解流體方程的過程中采用7個變量來表征物體幾何形狀，包括6個面積分數(shù)和單元的體積分數(shù)。面積分數(shù)為單元內(nèi)空白面積與總面積的比值，體積分數(shù)為單元內(nèi)空白體積與總體積的比值。同時，其能影響包括質(zhì)量、動量、能量守恒方程以及描述密度、壓力、溫度耦合情況的狀態(tài)方程在內(nèi)的所有控制方程形式，以此獲得準確的流體流場變化情況和自由液面運動情況。

1.3 Tru-VOF界面追蹤

自從1975年Hirt和Nichols提出Volume of Fluid(VOF)界面追蹤方法[12]以來，大多數(shù)的流體動力學(xué)計算軟件在追蹤自由界面形態(tài)變化時都會采用這一方法。研究所采用的Tru-VOF界面追蹤方法是在VOF方法基礎(chǔ)上進一步優(yōu)化而來的。該方法采用的是對流流體算法，不同于VOF方法計算空氣-流體的混合動力學(xué)性能，而是使用壓力或速度邊界代替氣體部分，不對氣體控制方程進行求解，準確追蹤流體運動過程中可能出現(xiàn)的尖銳界面，這樣既有效節(jié)約了計算時間，又一定程度地增加了模擬結(jié)果的精度。

Tru-VOF界面追蹤方法的主要控制方程如下：

(5)

式中，f是流體分數(shù)，其數(shù)值范圍是0～1之間，f=1代表著單元內(nèi)充滿流體，f=0表明單元內(nèi)為空氣；vf是擴散系數(shù)；FSOR是源項。

董爺爺不僅心靈手巧，還是學(xué)校里有名的環(huán)保達人。有一次冬天我吃完午飯，正巧碰到董爺爺在食堂的水槽邊用冰水在清洗碗筷，我就和他打招呼：“董爺爺好，食堂里不是有提供一次性的筷子嗎，您為什么不使用呢？”他笑著對我說：“用一次性的東西多不環(huán)保呀，而且，我自己帶碗筷也方便?！?/p>

1.4 泥沙模型

水流中泥沙的輸運方式主要分為推移質(zhì)運動和懸移質(zhì)運動，其中推移質(zhì)運動分為躍移質(zhì)和層移質(zhì)運動。泥沙模型中模擬泥沙顆粒輸運的示意圖如圖2所示。研究中泥沙沖刷過程的主要計算思路是基于沉積物輸運的基本理論，即當沙床剪切應(yīng)力τ高于臨界沙床剪切應(yīng)力τc時，沉積物沙床上的水流流動會推動沉積物顆粒開始運動。后續(xù)通過計算懸移質(zhì)、重力引起的泥沙沉降以及推移質(zhì)等來模擬沉積物沙床隨水流運動而產(chǎn)生的形態(tài)變化。

圖2 泥沙輸運方式示意圖

在水-沙交界面處，由于沙床剪切應(yīng)力以及漩渦的復(fù)合作用，泥沙往往會先被抬升，懸浮后沿水流方向輸運。鑒于難以通過流動動力學(xué)計算每一粒泥沙的運動軌跡，這里使用Mastbergen[13]等提出的經(jīng)驗?zāi)Ｐ汀ＥR界Shields數(shù)θcr,i的大小通過Soulsby-Whitehouse方程[14]計算：

(6)

沙床表面的粗糙度kS與泥沙的中值粒徑d50,packed成正比，

kS=croughd50,packed,

(7)

式中，crough是無量綱系數(shù)，通常為2.5。

懸移質(zhì)輸運速度的計算是基于經(jīng)驗公式：

(8)

其中，αi是泥沙被水流挾帶的系數(shù)，據(jù)經(jīng)驗?zāi)Ｐ腿?.018[13]，θi是局部希爾茲數(shù)，根據(jù)局部沙床剪應(yīng)力計算；ns是垂直于沙床表面向外的法向量；ulift,i是推移質(zhì)轉(zhuǎn)化為懸移質(zhì)的泥沙量。

在泥沙輸運的過程中，泥沙顆?？赡軙捎谧陨碇亓膽乙瀑|(zhì)中沉降到沙床上或在推移質(zhì)輸運的過程中靜止，從而產(chǎn)生泥沙淤積的現(xiàn)象。研究選取Soulsby[14]提出的沉降速度方程來計算泥沙沉降：

(9)

式中，υf是流體的運動粘度。

推移質(zhì)輸運的方式多為泥沙顆粒在水流沖刷作用下沿沙床滾動、跳躍，并在此過程中經(jīng)常與沙床表面泥沙進行交換。沙床單寬推移質(zhì)體積輸運率qb,i的計算采用Meyer-Peter[15]提出的公式：

(10)

式中，Φi為無量綱化的推移質(zhì)輸沙率；Φi=βi(θi-θcr,i1.5)，βi是無量綱系數(shù)，據(jù)經(jīng)驗?zāi)Ｐ腿?.0。

2 海底管道底部沖刷坑的形成和發(fā)展過程分析

為驗證模型本身的精確性與可靠性，選取Mao[1]等所做管道沖刷實驗為參照依據(jù)，簡化的管道沖刷數(shù)值水槽示意圖如圖3所示。水槽長6 m，高0.9 m，由于目的是研究二維管道沖刷過程，所以寬度取一個網(wǎng)格大小，即0.01 m。管道的材料參數(shù)設(shè)定為固壁(忽略其在沖刷過程中可能出現(xiàn)的形狀變化)，直徑D=0.1 m，在管道前方設(shè)置長40 D的泥沙底面，以使管道前方流場充分發(fā)展。為了保證出流不對沖刷過程產(chǎn)生影響，管道后方設(shè)置長20 D的泥沙底面。由于考慮海底滲流等因素對沖刷過程的影響，故將管道放置于沙床上方0.05 D。靜水深0.35 m，水流流向從左至右，進口的邊界條件設(shè)置為速度邊界，水流初始速度為0.25 m/s。同時，為了不讓泥沙在入水口處就直接受到水流沖刷作用，在上游入口設(shè)置了厚0.01 m，高0.3 m的穩(wěn)水擋板。下游出口邊界條件選擇壓力邊界以設(shè)定流體高度保持模擬過程中水位保持不變；底部采用無滑移壁面條件，粗糙度為2.5 d50；頂部選取壓力邊界，模擬自然環(huán)境中大氣壓對水流流動的影響；其余均為對稱邊界。數(shù)值水槽的底部鋪設(shè)厚0.3 m的沙床，泥沙的密度為2 700 kg/m3，中值粒徑為0.19 cm，Shields數(shù)的臨界值是0.048，孔隙度為0.4，休止角為32°。數(shù)值水槽的整體網(wǎng)格數(shù)為116 880，網(wǎng)格大小為0.01 m。為了獲得管道沖刷過程中更加精確的數(shù)據(jù)，在管道前后特別加密了網(wǎng)格(圖3虛線之間部分)，網(wǎng)格大小為0.005 m。

圖3 管道沖刷數(shù)值水槽示意圖

將數(shù)值計算結(jié)果與Mao[1]等的實驗結(jié)果進行了比較，沙床高程結(jié)果的對比、模擬過程中流場的壓強及速度分量變化如圖4所示。圖4a、圖4b是研究沙床高程數(shù)值計算與Mao[1]等實驗10 min和30 min的結(jié)果對比圖。從兩者曲線整體來看，數(shù)值計算與實驗結(jié)果呈現(xiàn)出較好的一致性。兩者之間存在些許差異是由于管道并未直接放置于沙床表面，而是存在一定間隙，所以管道后方堆積體位置較實驗結(jié)果略微偏右，但隨時間推移該間隙對整體堆積體的位置變化影響可以忽略。同時，由于未考慮海底滲流效應(yīng)對泥沙顆粒的浮力作用、管道系統(tǒng)以及在沖刷過程中管道表面形狀變化對沖刷過程的影響，故沖刷坑深度以及堆積體的高度會比實驗結(jié)果略高。以上說明，研究數(shù)值模型對于海底管道底部沖刷問題具有較好的適用性和較高的計算精度。

圖4c、圖4d分別為沖刷時間10 min和30 min的流場流線圖。從圖4c、圖4d中可以看出，水流在到達管道之前為均勻來流，靠近沙床的水流速度較低。在遇到管道時，水流出現(xiàn)分流，在管道上下兩側(cè)水流速度達到最大，最大值為0.406 m/s。這幾乎是入流速度的1.6倍，導(dǎo)致了在管道底部發(fā)生射流現(xiàn)象，并因此產(chǎn)生沖刷。在管道后上方，水流速度近似或等于0，形成渦流區(qū)，部分懸移質(zhì)沉降，推移質(zhì)靜止，導(dǎo)致在管道后方逐漸形成淤積。而當淤積泥沙體積逐漸變大時，管道后下方水流速度緩慢變小，堆積體右側(cè)水流速度為負，形成漩渦，減緩堆積泥沙輸運。在管道的下游側(cè)，沙床高程與初始高程相近，沒有觀察到?jīng)_刷現(xiàn)象。這是因為沙床上坡度較大的沖刷坑會使水流幾乎垂直地從管道下方?jīng)_刷而出，避免了漩渦與沙床的相互作用。

圖4e、圖4f是模擬時間10 min和30 min時的流場壓強圖。由圖4e、圖4f中可以看出，由于穩(wěn)定流作用，管道底部的壓強始終較大，為3.6 kPa左右。這也解釋了為何在沖刷過程中管道底部沖刷坑深度逐漸增加。結(jié)合圖4c、圖4d以及10 min到30 min的壓強均為正壓可以看出，相同條件下，隨著沖刷時間的增加，管道底部沖刷坑加深，水流流速減小，管道底部沙床所受壓強增大。該現(xiàn)象也符合伯努利原理，從側(cè)面驗證了數(shù)值模型的可靠性。

圖4g、圖4h是10 min和30 min時的z方向速度分量圖。由圖4g、圖4h中可以看出，在遇到管道之前，水流z方向速度分量為0。遇到管道時，水流分流，部分向上速度達到z方向最大值，部分向下沖刷管道底部沙床，形成沖刷坑，并開始沿水流方向輸運泥沙顆粒。雖然部分水流沿沖刷坑坡度流出會給沙床表面泥沙一個z方向的速度，但由于管道后方水流合流，該速度很快為0，導(dǎo)致懸浮泥沙顆粒逐漸沉降形成淤積。這也是管道后方會形成泥沙堆積體的原因之一。

圖4 管道沖刷數(shù)值結(jié)果分析圖

從以上結(jié)果分析圖及驗證部分中可以清楚得出穩(wěn)定流作用下海底管道底部沖刷坑形成和發(fā)展的原因。水流在到達管道之前，流場穩(wěn)定為均勻流，流速大小為設(shè)定值。水流在遇到管道時分流，在管道上下兩側(cè)水流速度達到最大，并在管道底部出現(xiàn)射流現(xiàn)象，產(chǎn)生沖刷坑。射流沿沖刷坑坡度流出會給沙床表面泥沙一個z方向的速度，但由于管道后方合流，水流產(chǎn)生速度差形成漩渦，導(dǎo)致部分懸移質(zhì)沉降，推移質(zhì)靜止形成淤積。同時，因為沙床表面沖刷坑坡度較大，射流幾乎垂直從管道下方?jīng)_出，防止了管道后方由于速度差形成漩渦與沙床再次相互影響，所以只有管道附近沙床會受到影響。

3 海底管道底部沖刷的影響因素

為了更深入地研究海底管道底部的沖刷過程，通過調(diào)整驗證模型中的流速、靜水深度、管道直徑以及管道與沙床間隙等參數(shù)，進一步探討相關(guān)參數(shù)對沖刷坑發(fā)展過程產(chǎn)生的影響。

3.1 流速

圖5 不同入流流速條件下數(shù)值結(jié)果

3.2 水深

在其他參數(shù)設(shè)置與驗證算例相同時，設(shè)置不同的靜水水深，海底管道底部沖刷過程的數(shù)值結(jié)果如圖6所示。圖6a、圖6b分別是模擬時間為10 min、30 min時，在靜水深為0.2 m、0.3 m、0.4 m條件下的沙床高程數(shù)值結(jié)果圖。當水深由0.2 m到0.4 m時，管道底部沖刷坑的深度逐漸減小，管道后方堆積體體積縮小，離管道的距離逐漸變小。圖6c、圖6d為30 min時，水深為0.2 m和0.4 m的流場z方向速度分量圖。從圖6c、圖6d可以看出，隨著靜水深度的增加，水流在經(jīng)過管道截面時的流速減小，管道對周圍流場的影響逐步減弱。水流分流后，部分向下方?jīng)_刷沙床的水流z方向速度分量減小，管道后方射流z方向的流速減小以致水流攜帶泥沙的能力降低，沙床表面沖刷坑深度變淺。同時，水流經(jīng)管道后合流速度加快，管道后方淤積泥沙輸運速率減緩，堆積體離管道的距離更近。這也就是圖6a、圖6b兩幅圖中沙床高程變化的原因。

圖6 不同靜水深度條件下管道沖刷過程的數(shù)值結(jié)果

3.3 管道直徑

在其他參數(shù)設(shè)置與驗證算例相同時，設(shè)置不同的管道直徑，海底管道底部沖刷過程的數(shù)值結(jié)果如圖7所示。圖7a、圖7b分別為模擬時間是10 min、30 min時，在管道直徑是0.05 m、0.1 m、0.2 m條件下的沙床高程數(shù)值結(jié)果圖。當管徑為0.05 m(即靜水深與管道直徑的比值h/D=7)時，無量綱化(S/D即沖刷深度/管道直徑)后的沖刷深度會比其余兩種情況略深，堆積體也會更向右移。而當管徑為0.1 m、0.2 m(即h/D≤3.5)時，沙床高程無量綱化后的結(jié)果極為相近，沖刷坑深度、范圍及形狀相似，且隨數(shù)值模擬時間增長，曲線更為接近。圖7c、圖7d為30 min時，管道直徑是0.05 m和0.2 m的管道沖刷流場速度圖。從圖7c、圖7d可以看出，雖然相同水深條件下，隨著管徑的增大，水流經(jīng)過管道橫截面時前后流場所受到的影響更大，管道下方水流流速增大，沖刷坑深度也會更深，但是由于管徑不同，將沖刷深度無量綱化(S/D)后則可以看出，管徑小的沖刷深度反而會更大。

圖7 不同管道直徑條件下管道沖刷過程的數(shù)值結(jié)果

3.4 管道與沙床間隙

在其他參數(shù)設(shè)置與驗證算例相同時，設(shè)置不同管道與沙床的間隙值，海底管道底部沖刷過程的數(shù)值結(jié)果如圖8所示。圖8a、圖8b分別為模擬時間是10 min、30 min時，在管道與沙床間隙是0.005 m、0.02 m、0.04 m條件下的沙床高程數(shù)值結(jié)果圖。由圖8可知，隨著管道與沙床間隙的增大，沙床沖刷坑的深度在逐漸減小，管道沖刷的階段可以快速從射流期過渡到尾流期[2]。從間隙0.005 m到間隙0.02 m的曲線對比中可以看出，隨著管道上移，水流在經(jīng)過管道時流場的變化也會隨之上升，導(dǎo)致沖刷深度減小，產(chǎn)生的射流現(xiàn)象減弱，泥沙輸運能力降低以致堆積體積逐漸增大。但從間隙0.02 m到間隙0.04 m的沙床高程曲線對比中可以看出，沖刷坑深度變小，但由于管道與沙床間隙過大，水流流速較快，使得堆積體表面泥沙快速輸運，以致間隙0.04 m時的堆積體積較間隙0.02 m的更小。圖8c、圖8d是30 min時，管道與沙床間隙是0.005 m和0.04 m的流場z方向速度分量圖。從圖8c、圖8d中可以看出，隨著管道的上移，管道下方?jīng)_刷沙床水流的z方向速度則會降低，輸運泥沙能力減弱，導(dǎo)致沖刷坑深度逐漸減小。

圖8 不同管道與沙床間隙條件下管道沖刷過程的數(shù)值結(jié)果

4 結(jié)論