陳海登，陶祥海，李 玲

(廣東電網(wǎng)有限責(zé)任公司湛江供電局，廣東 湛江 524005)

海底管線在海洋工程中得到了普遍應(yīng)用。當(dāng)管線鋪設(shè)在沙質(zhì)海床上時，管線的存在改變了周圍的流場，增加了周圍沉積物的輸送速率，并直接導(dǎo)致管線周圍的局部沖刷，從而引起海洋工程結(jié)構(gòu)的破壞。目前許多工作致力于局部沖刷數(shù)值模型的開發(fā)，采用不同的湍流模型結(jié)合不同的經(jīng)驗(yàn)輸沙公式模擬沖刷進(jìn)程，并用于預(yù)測管線下方的最大沖刷深度、管線上下游海床變化和漩渦脫落情況。Leeuwestein等[1]采用k-ε模型和推移質(zhì)輸移經(jīng)驗(yàn)公式模擬海床變化，但是模型中沒有考慮懸移質(zhì)的輸運(yùn)，導(dǎo)致沖刷坑計算值偏大。van Beek等[2]基于k-ε湍流模型，建立了一個只考慮懸移質(zhì)輸運(yùn)的沖刷模型，并將該模型用于預(yù)測管線下方的沖刷，預(yù)測的沖刷坑與試驗(yàn)結(jié)果吻合較好，但是下游沖刷變化有所偏差。Br?rs[3]建立了一個既考慮密度影響又考慮懸移質(zhì)和推移質(zhì)輸運(yùn)的沖刷模型，利用該模型預(yù)測了清水沖刷，結(jié)果與Mao[4]的試驗(yàn)結(jié)果吻合較好，但是該模型沒有模擬出管線后方周期性的漩渦脫落。Li等[5]采用大渦模擬方法建立了管線周圍局部沖刷的數(shù)值模型，模型中假定河床剪切應(yīng)力與清水沖刷下泥沙臨界剪切應(yīng)力相等，或直接采用動床沖刷下不受擾流的床面剪切應(yīng)力，預(yù)測的最大沖刷深度和沖刷坑形狀與試驗(yàn)數(shù)據(jù)吻合良好。該模型的主要優(yōu)點(diǎn)是不需要任何泥沙輸移公式，但是不能準(zhǔn)確地預(yù)測沖刷過程。Liang等[6-7]建立了一個數(shù)值模型，用于預(yù)測清水和動床條件下管線周圍局部沖刷的時間發(fā)展歷程。該模型同時考慮了懸移質(zhì)和推移質(zhì)的輸運(yùn)，在曲線坐標(biāo)系下采用有限差分法進(jìn)行求解。在模型中使用了新的時間推進(jìn)算法和沙滑模型，分別采用標(biāo)準(zhǔn)k-ε模型和大渦模擬計算了整個沖刷過程。結(jié)果表明，標(biāo)準(zhǔn)k-ε模型的計算結(jié)果比大渦模擬計算結(jié)果更加精確，沖刷剖面對流場的計算網(wǎng)格密度的依賴性較小。

本文通過有限元法[8]建立了一個沖刷數(shù)值模型，基于非定常雷諾平均N-S方程和標(biāo)準(zhǔn)k-ε湍流模型來模擬水流流動，采用經(jīng)驗(yàn)公式來模擬推移質(zhì)運(yùn)動，利用動網(wǎng)格方法[9]來捕捉水沙交界面，所有方程采用特征線法[10]進(jìn)行時間離散。通過高雷諾數(shù)下二維鈍體繞流和Mao[4]的清水沖刷試驗(yàn)對模型的有效性進(jìn)行驗(yàn)證。

1 數(shù)值模型

1.1 水流控制方程

流場求解采用非定常雷諾平均N-S方程，湍流模型選擇k-ε模型，控制方程如下：

(1)

(2)

(3)

(4)

式中：ui、uj分別為在xi、xj(i=1,2,3;j=1,2,3)方向上的水流速度；ρ為水的密度；p為壓力；t為時間；ν為動力黏度；τij為雷諾應(yīng)力；k為湍動能；νt為湍流黏度，νt=Cμk2/ε；ε為湍動能耗散率；δij為克羅奈克函數(shù)(當(dāng)i=j時，δij=1;當(dāng)i≠j時，δij=0)；Cμ為經(jīng)驗(yàn)常數(shù)，取值為0.09；Pk為由于平均速度梯度引起的湍動能的產(chǎn)生項；σk、σε、C1ε、C2ε均為經(jīng)驗(yàn)常數(shù)，分別取值為1.0、1.3、1.44、1.92。

為了精確模擬沖刷，在求解過程中使用了動網(wǎng)格方法。動網(wǎng)格方法的基本思想是在每個時間步內(nèi)，通過對流體區(qū)域的網(wǎng)格進(jìn)行更新以實(shí)現(xiàn)由于邊界運(yùn)動引起的計算域與其物理量的變化，采用任意拉格朗日-歐拉(ALE)方法描述流體運(yùn)動的控制方程，通過插值運(yùn)算從舊網(wǎng)格映射得到網(wǎng)格更新后的各個物理量，然后對每個時間步內(nèi)網(wǎng)格上的物理量進(jìn)行迭代求解來獲得流動的動態(tài)演化結(jié)果[11]。當(dāng)前主要的動網(wǎng)格方法有虛擬結(jié)構(gòu)法、偏微分方程法、代數(shù)法、網(wǎng)格重構(gòu)法和網(wǎng)格變形與局部或全場網(wǎng)格重構(gòu)結(jié)合[12]。本文通過構(gòu)建方程，求解偏微分方程計算網(wǎng)格節(jié)點(diǎn)坐標(biāo)變化。

采用動網(wǎng)格方法后， N-S方程中連續(xù)方程保持不變，僅需將動網(wǎng)格速度添加進(jìn)動量方程、湍動能、湍動能耗散率方程中：

(5)

(6)

(7)

式中：vj為xj(j=1,2,3)方向上動網(wǎng)格速度。

1.2 沖刷模型

導(dǎo)致沖刷的主要原因是懸移質(zhì)和推移質(zhì)的運(yùn)動[13]。懸移質(zhì)輸移模型適合用于挾沙運(yùn)移，而不適用于只有沉積物形成的海床。本文只考慮推移質(zhì)運(yùn)移導(dǎo)致的沖刷。由于van Rijn推移質(zhì)輸沙率公式[14]中各參數(shù)意義明確，計算方便，易于理解，可以直接運(yùn)用于河流、河口和海岸等地區(qū)，且本文分別采用有限元法和特征線法進(jìn)行空間和時間離散，可以避免不穩(wěn)定問題[15]，本文采用van Rijn公式計算推移質(zhì)輸沙率：

qb=ψu(yù)b

(8)

τωcr0=(ρs-ρ)gd50θcr

β=tan-1‖Hzb‖

式中：qb為沿著海床單寬推移質(zhì)通量；ψ為海床上的沉積物濃度；ub為平行于海床上推移質(zhì)層面的水流流速；δb為推移質(zhì)層的厚度，取δb=2d50；d50為顆粒的中值粒徑；ρs為砂子的密度；T為海床面無量綱剪切力；τω和τωcr分別為海床面剪切力和臨界剪切力，D*為無量綱粒徑；g為當(dāng)?shù)刂亓铀俣?；τωcr0為平整床面臨界剪切力；φ為顆粒的休止角；β為顆粒最大的滑動角；α為推移質(zhì)層面流速ub與最大滑動坡度方向的夾角；θcr為水平面上的Shields常數(shù)，一般取0.048；H為水平梯度算子；zb為床面演變高度，采用Exner方程計算；n為海床的孔隙率。

1.3 離散方法

采用有限元法對該模型進(jìn)行空間離散，特征線法對時間進(jìn)行離散。一般地，任意守恒型方程可以寫作如下形式：

(9)

式中：Φ為未知量；F為對流通量；G為擴(kuò)散通量；Q為源項。式(1)(5)(6)(7)都可寫成式(9)形式。

除去三階項，將式(9)沿特征線方向展開，采用張量記法的形式，m時刻的變化量為

ΔΦ=-Δt(·F+·G+Q)m+

(10)

式中：Δt為時間步長；u為速度項；ΔΦ為未知量Φ變化量，同時，ΔΦ滿足

(11)

(12)

式中：Ω為計算域；n為界面外法線；Γ為積分域邊界。

2 模型驗(yàn)證

2.1 雷諾平均的非線性k-ε湍流模型驗(yàn)證

采用笛卡爾坐標(biāo)系下二維方形鈍體繞流對雷諾平均的非線性k-ε模型進(jìn)行驗(yàn)證，模型總長度為18H(H為方形鈍體寬度)，高度為2H，方形鈍體位于模型左端3H處，模型參數(shù)及邊界情況見文獻(xiàn)[16]。模型左側(cè)為流體速度進(jìn)口，右側(cè)為自由出流，上下設(shè)置為滑移邊界。將模擬結(jié)果與文獻(xiàn)[16-18]結(jié)果進(jìn)行了對比分析。設(shè)定雷諾數(shù)Re=40 000，與文獻(xiàn)[17]中一致。

圖1為鈍體左側(cè)0.2H處的水平速度和豎向速度分布與文獻(xiàn)[16-18]結(jié)果的對比，其中u0為入口流速。由圖1可知，計算結(jié)果與文獻(xiàn)結(jié)果較為吻合，變化趨勢一致，而且在近壁區(qū)域本文模擬結(jié)果更接近于實(shí)測結(jié)果，精度更高。

圖1 速度分布對比

圖2為計算得到的流線圖，從圖中可以看出，本文建立的模型可以捕捉到鈍體背流面的主回流區(qū)、頂部的分離區(qū)、迎流面的小渦旋結(jié)構(gòu)及背流面角點(diǎn)處的小渦旋。相比于文獻(xiàn)[16]，本文模擬的流場可以更好地展現(xiàn)出背流面角點(diǎn)處的小渦旋，精度更高。本文計算出的主回流區(qū)長度介于鈍體后方(1.0H～7.5H)范圍內(nèi)，與文獻(xiàn)[18]的實(shí)測結(jié)果以及文獻(xiàn)[16-17]的模擬結(jié)果相一致，表明本文計算方法具有較高的準(zhǔn)確性，可以用于沖刷模型的驗(yàn)證。

圖2 流場結(jié)構(gòu)

2.2 沖刷模型驗(yàn)證

為了驗(yàn)證模擬結(jié)果，采用Mao[4]的試驗(yàn)數(shù)據(jù)進(jìn)行對比分析，其試驗(yàn)中管線沖刷條件為未擾動的希爾茲常數(shù)θ=0.048。在Mao[4]的試驗(yàn)中，管徑D=100 mm，水深h=3.5D。海床表面顆粒的中值粒徑d50=0.36 mm，進(jìn)口水流速度為350 mm/s。

圖3 網(wǎng)格劃分

根據(jù)此試驗(yàn)，建立了二維有限元沖刷模型，上游計算長度為10D，下游計算長度為20D，水深為3.5D。為了避免在計算過程中網(wǎng)格拓?fù)浣Y(jié)構(gòu)的變化，在管線下方設(shè)置了一個初始沖刷坑，深度為0.1D，網(wǎng)格劃分如圖3所示。邊界條件設(shè)置描述如下：① 進(jìn)口(圖3左側(cè))：y向速度設(shè)置為0，x向給定水平速度u1、湍動能k和耗散率ε，這些值是根據(jù)事先建立的一個相同水深，長度為100D的明渠模型計算獲得，將此明渠的出口作為本模型的進(jìn)口。②出口(圖3右側(cè))：邊界設(shè)置為自由出流，即速度、湍動能和耗散率的梯度為0，并給定壓力參考點(diǎn)。③模型頂部設(shè)置為剛蓋，即速度等變量在y向的梯度均為0。④在海床和管線表面采用了壁面函數(shù)[19]。

根據(jù)模擬結(jié)果，討論在初始沖刷階段(t≤30 min)的河床變化，并與Mao[4]的清水沖刷試驗(yàn)結(jié)果進(jìn)行了對比。圖4顯示了不同時刻的沖刷形態(tài)。如圖4(a)所示，開始階段，由于管線上下游形成的壓差，海床表面出現(xiàn)滲流。在滲流出口處，當(dāng)壓力超過泥沙的浮重時，沖刷開始發(fā)生。隨著滲流的進(jìn)一步發(fā)展，大量泥沙向下游沖去，沉積在管線后方，如圖4(b)所示。經(jīng)過一段時間后，管線和海床之間的間隙變大，并且泥沙在管線后方形成沙脊，如圖4(c)所示。此過程持續(xù)，沖刷深度不斷增加，沙脊也向下游移動，如圖4(d)所示。

圖4 不同時刻的海床沖刷形態(tài)

圖5 不同時刻局部沖刷對比

圖5為本文模擬結(jié)果與Mao[4]的試驗(yàn)數(shù)據(jù)和Liang等[7]的數(shù)值模擬結(jié)果的對比，可以得出：

a. 本文數(shù)值模型的計算結(jié)果與Mao[4]試驗(yàn)數(shù)據(jù)非常一致，尤其是在管線下方與沙脊之間，本文模擬的沖刷形態(tài)比Liang等[7]的k-ε模型和大渦模擬計算的結(jié)果更準(zhǔn)確。

b. 從堆積形態(tài)上看，Liang 等[7]模擬的堆積高度高于Mao[4]的試驗(yàn)結(jié)果，而本文模擬值略低于試驗(yàn)結(jié)果，但更接近于試驗(yàn)結(jié)果。

c. 在沙脊后部，模擬的結(jié)果與試驗(yàn)相比有些不同，可能是由于未考慮泥沙休止角和忽略懸移質(zhì)輸運(yùn)。如圖5(a)所示，在沖刷初始階段，模擬結(jié)果與試驗(yàn)結(jié)果吻合非常好，原因是此時的懸移質(zhì)濃度很低，沖刷發(fā)生主要以海床推移質(zhì)為主。隨著時間的推移，懸移質(zhì)的濃度開始增加，沖刷變成由懸移質(zhì)和推移質(zhì)共同作用，這就導(dǎo)致了模擬結(jié)果與試驗(yàn)結(jié)果有所差異(圖5(b))。沙脊的模擬結(jié)果與Liang等[7]的大渦模擬結(jié)果相似，因?yàn)槟Ｐ投紱]有考慮休止角，但是在起始階段，Liang等[7]的大渦模擬高估了沖刷深度(圖5(a))。

3 結(jié) 語

本文建立了一個有限元數(shù)值模型，用于預(yù)測海底管線在穩(wěn)定水流下的局部沖刷。采用標(biāo)準(zhǔn)k-ε湍流封閉模型求解非定常雷諾平均N-S方程來模擬流場，采用有限元對模型進(jìn)行空間離散，特征線法對時間進(jìn)行離散。模擬的二維方形鈍體繞流結(jié)果與文獻(xiàn)[16-18]的結(jié)果吻合較好，驗(yàn)證了湍流模型的有效性。模擬的沖刷坑形狀和沖刷深度與Mao[4]的試驗(yàn)結(jié)果吻合較好，在沖刷初始階段的模擬結(jié)果更接近于試驗(yàn)值。與文獻(xiàn)[7]的模擬結(jié)果相比，管線下方與沙脊之間的海床演變得到了更精確的模擬，但是沙脊形成后的模擬結(jié)果與試驗(yàn)結(jié)果有所差異，可能是由于模擬中未考慮休止角和懸移質(zhì)的作用。

致謝：本文相關(guān)研究工作得到河海大學(xué)水利水電學(xué)院趙蘭浩教授和廣東省電力設(shè)計研究院李健華博士的幫助，特此致謝。