海床破壞下的海底管道路由區(qū)水動力沖淤特性研究與應(yīng)用

淳明浩，于德周，楊肖迪，羅小橋，姚志廣，徐 爽

1.中國石油集團工程技術(shù)研究有限公司，天津 300451

2.中國石油集團海洋工程重點實驗室，天津 300451

3.中國石油集團海洋工程（青島）有限公司，山東青島 266520

海底管道作為海底油氣資源集輸?shù)闹匾獦屑~，其安全運維是海上油氣資源開發(fā)中的關(guān)鍵環(huán)節(jié)之一[1-4]。海底管道路由區(qū)在未受到外來破壞之前，通常會處于一種相對的海床沖淤平衡狀態(tài)[5-7]。然而，由于管道所處路由區(qū)的水動力條件變化、人為作業(yè)活動或自然因素破壞導(dǎo)致的周邊海洋環(huán)境條件改變，造成原有沖淤平衡狀態(tài)的破壞[5-6]，均可能引起海底管道路由區(qū)局部或整體的沖刷侵蝕[3-4]，導(dǎo)致海底管道上覆土體保護層減薄和侵蝕[4-6]，威脅海底管道安全。

管道路由區(qū)海床的穩(wěn)定是海底管道安全與穩(wěn)定的有力保障，諸多海底管道安全風(fēng)險常由海床的沖淤變化引起[5-7]。目前，針對海底管道及其路由區(qū)的研究多集中于管道本體、路由區(qū)選址等方面，對于管道路由區(qū)海床沖淤侵蝕的研究較少[1-6]。

同時，海底管道路由區(qū)沖淤測量與分析的方法主要來源于設(shè)備觀測、原型測量分析、物理模型及數(shù)值模型模擬等手段[8-11]。本文通過研究，提出基于海底管道路由區(qū)水動力沖淤特征的分析方法，對于明確水動力條件變化對管道路由區(qū)的影響可以起到支撐與促進作用，可與海底管道路由區(qū)的觀測調(diào)查形成綜合互補的路由區(qū)風(fēng)險分析評估技術(shù)，有助于提高海底管道路由區(qū)相關(guān)風(fēng)險的分析評估效率和管道運維保護的科學(xué)性。

1 研究區(qū)域概況

某海底管道位于廣東省深圳市境內(nèi)珠江口海域，起點為深圳大鏟島，終點為深港邊界，管道途經(jīng)采砂區(qū)、銅鼓航道、礬石淺灘區(qū)、臨近錨地等特殊地段，海底段全長14.78 km[1-2]（以起點為KP 0、終點為KP 14.78表示）。2018年12月探測調(diào)查發(fā)現(xiàn)，管道路由區(qū)受到海底采砂作業(yè)破壞，形成巨大采空區(qū)，對原有水動力沖淤平衡造成了破壞[2]。為研究海底管道路由區(qū)當(dāng)前水動力沖淤特性，在路由區(qū)與采砂區(qū)勘測詳查的基礎(chǔ)上，2019年進一步開展了水動力參數(shù)觀測、表層沉積物取樣與組成測試，根據(jù)獲得的水動力參數(shù)和沉積物組成參數(shù)，研究建立水動力沖淤數(shù)值模擬模型，評價海底管道路由區(qū)的沖刷侵蝕趨勢，為管道運行維護和風(fēng)險評估提供支持。海底管道地理位置[1]如圖1所示，水動力參數(shù)觀測與沉積物取樣位置分布[2]如圖2所示。

圖1 海底管道地理位置

圖2 水動力參數(shù)觀測與海底沉積物取樣位置分布示意

2 水動力沖淤特性方法構(gòu)建

2.1 原理與技術(shù)路線

本文所研究管道路由區(qū)海床由于采砂破壞而形成了巨大凹坑，因而其原水動力沖淤平衡已破壞，當(dāng)前水體的紊流狀態(tài)對海床沖淤變化影響最大[7-8]，因此采用紊流模型進行沖淤分析。

在水動力模擬分析中，當(dāng)采用紊流模型模擬流場時，可獲得任意水質(zhì)點的流速，計算海床泥沙受到的剪切力，再根據(jù)土質(zhì)的泥沙起動剪切力，判斷海床表面的泥沙是否運移起動[12]。模擬分析過程中，通常將運動的泥沙分為懸移質(zhì)和推移質(zhì)，由海床表面的剪切力和推移質(zhì)輸運公式得到推移質(zhì)的輸沙率，通過求解懸移質(zhì)質(zhì)量濃度控制方程得到分析區(qū)域的泥沙質(zhì)量濃度，進而計算懸移質(zhì)的輸沙率，總的輸沙率為推移質(zhì)輸沙率和懸移質(zhì)輸沙率之和[8-9]。在包括海床面的封閉空間內(nèi)，利用質(zhì)量守恒原理求得海床面的變化，經(jīng)過多次的迭代，得到海床穩(wěn)定后的形態(tài)[14-15]。針對本研究區(qū)域的管道路由區(qū)水動力沖淤模擬，共建立了流場模型、泥沙輸運模型及地形變化模型三種模型[10-12]。本次水動力沖淤模擬模型的基本邏輯關(guān)系如圖3所示。

圖3 海底水動力沖淤模型邏輯關(guān)系

2.2 數(shù)值模型研究

根據(jù)上述確定的紊流水動力沖淤模型，確定本次水動力沖淤模型的各個方程組成。

其中，模擬管道路由區(qū)的流場模型由時均雷諾方程和連續(xù)方程組成[7-9]。在笛卡爾坐標(biāo)系下的無因次時均雷諾方程如下[14-17]：

式中：t為時間，s；x、y、z是笛卡爾坐標(biāo)系坐標(biāo)，m；η為水面標(biāo)高，m；p為水體壓強，Pa；u、v、w分別對應(yīng)x、y、z坐標(biāo)方向上的水體速度分量，m/s；Re為雷諾系數(shù)；τ11為xx平面的剪切應(yīng)力，Pa；τ12為xy平面的剪切應(yīng)力，Pa；τ22為yy平面的剪切應(yīng)力，Pa。

無因次時均連續(xù)方程[14-17]如下：

由于以上基本方程構(gòu)成的方程組并不封閉，本次數(shù)值模擬采用標(biāo)準(zhǔn)k-ε紊流方程進行封閉。標(biāo)準(zhǔn)k-ε紊流方程[14-17]如下：

表1 標(biāo)準(zhǔn)k-ε紊流模型的參數(shù)取值

表1 中的Cμ、σk、σε、Cε1、Cε2均為紊流模型經(jīng)驗常數(shù)。

2.3 泥沙輸運模型研究

泥沙運動方式通常分為懸移質(zhì)移動和推移質(zhì)移動兩種[14,15,19]。在水動力模擬分析中，運動狀態(tài)泥沙中的推移質(zhì)為緊貼底面很薄的一層泥沙顆粒，其在運動中不離開海床底面，而懸移質(zhì)會因?qū)α骱蛿U散而分布到整個流場中[17-19]。根據(jù)2019 年在如圖2所示的管道路由區(qū)進行的海底沉積物取樣成分分析報告，本研究區(qū)的沉積物組成類型為以淤泥質(zhì)砂質(zhì)黏土和砂質(zhì)黏土為主，還有較少粉砂，故本次研究中對于泥沙運動同時考慮泥沙懸移質(zhì)和推移質(zhì)的運動。

2.3.1 懸移質(zhì)輸運模型分析

本研究區(qū)懸移質(zhì)的輸運受泥沙濃度控制[20-22]。本次模擬模型采用對流擴散方程來求解泥沙質(zhì)量濃度，對流擴散方程表示如下[7-9,20]：

式中：c為泥沙的質(zhì)量濃度，g/m3；ωs為泥沙顆粒在水中的沉降速度，m/s；σc為紊流Schmidt數(shù)，表示泥沙擴散系數(shù)與渦黏系數(shù)的關(guān)系[19-20]，本研究中σc=1.0。

海床附近的泥沙質(zhì)量濃度求解公式[7-9,20]如下：

式中：d50為靜水深度，m。Δb為海床的參考高度，即海床到懸移質(zhì)輸沙的距離，m。u*為泥沙摩阻速度，N·s/m2；u*cr為泥沙的臨界摩阻速度，N·s/m2。D*為水體黏滯系數(shù)[20-21]。

本次模擬模型的泥沙質(zhì)量濃度在海底的邊界條件設(shè)定為：當(dāng)T＞0時，利用式（7）計算海底邊界的泥沙質(zhì)量濃度；當(dāng)T≤0時，海底泥沙的質(zhì)量濃度在邊界法線方向梯度為0。在運動水體的圓柱表面和自由面，水體邊界法線方向的泥沙流量為0，泥沙質(zhì)量濃度計算方程[7-9,20]如下：

式中：β為水體邊界與垂線之間的夾角，（°）；n為水體邊界法線方向的單位向量[22]。

求解方程（8），得到分析區(qū)域任一點的泥沙質(zhì)量濃度，結(jié)合流場分析結(jié)果，懸移質(zhì)輸沙率qs的計算公式[7-9,20]如下：

式中：yb為海床的垂直坐標(biāo)，m；H為b點的水體深度，m。

2.3.2 推移質(zhì)輸運模型分析

本研究區(qū)推移質(zhì)的輸運受泥沙密度、流體密度及流體速度控制[20-22]，本次數(shù)值模擬中推移質(zhì)輸沙率[23-26]的計算公式[7-9,20]如下：

式中：qb為單位寬度內(nèi)的體積輸沙率，kg/s；g為重力加速度，m/s2；s=ρs/ρ，為泥沙密度和流體密度的比值。

2.4 地形變化模型設(shè)定

本研究區(qū)的海底地形變化模型的建立依據(jù)為泥沙在海床區(qū)域的質(zhì)量守恒原則[25-28]，因此本次模型對于地形變化引起的海床變化采用的控制方程[7-9,20,29]如下：

式中：po為泥沙孔隙率，%；qT為總泥沙輸運率，kg/s，等于推移質(zhì)輸沙率qb和懸移質(zhì)輸沙率qs之和[18-22]。

3 數(shù)值模型分析計算

由于管道路由區(qū)受到人為作業(yè)的破壞，形成了海底采砂坑，造成管道路由區(qū)域海床形態(tài)的巨大變化，這種地形變化對管道路由區(qū)的水動力特征產(chǎn)生了影響，導(dǎo)致路由區(qū)不同區(qū)域的泥沙沖刷侵蝕狀況的改變。為深入分析這種沖淤變化及其影響，采用不同模型尺度進行數(shù)值模擬計算。

為了提高數(shù)值模擬計算的準(zhǔn)確度以及更好地模擬分析不同尺度下的路由區(qū)沖淤特性，對整個研究區(qū)設(shè)置了大區(qū)域模型和小區(qū)域模型。其中，大區(qū)域模型以涵蓋管道路由區(qū)、外部海域為基礎(chǔ)，主要模擬珠江口外部的海洋水動力對路由區(qū)的影響。小區(qū)域模型以管道路由區(qū)、周圍徑流區(qū)及潮流區(qū)為基礎(chǔ)，主要模擬潮流和徑流水動力對路由區(qū)的影響。

3.1 大區(qū)域模型設(shè)置

本次水動力沖淤模擬的大區(qū)域模型的計算范圍為針對海底管道所處的廣東省深圳市珠江口外的近岸海域，計算區(qū)域范圍為12.7°N～29.4°N、105.6°E～124.5°E，并對工程模型區(qū)域網(wǎng)格進行加密。該模型的網(wǎng)格布置與地形分布如圖4所示。

圖4 大區(qū)域模型計算范圍與網(wǎng)格劃分圖（黑色區(qū)域為研究區(qū)）

3.2 小區(qū)域模型設(shè)置

本次水動力沖淤模擬的小區(qū)域模型的計算范圍為海底管道路由區(qū)臨近的珠江口伶仃洋和外部海域。本研究中，小區(qū)域模型采用三維垂向平均模式。為了更好擬合水-陸邊界，模型利用非結(jié)構(gòu)三角網(wǎng)格技術(shù)對計算區(qū)域進行空間剖分。此類型三角網(wǎng)格設(shè)計靈活，便于調(diào)整網(wǎng)格疏密，有助于減少鋸齒狀岸線對模擬計算結(jié)果的不利影響。小區(qū)域模型計算范圍與地形分布如圖5所示。

圖5 小區(qū)域模型計算范圍與水深地形分布圖（圖中黑線為管道位置）

3.3 邊界條件設(shè)置

由于管道路由區(qū)域受到外海來水、內(nèi)陸來水、海底沉積物等多重影響，因此采用多重邊界條件對模型進行限定，包括開邊界條件、閉邊界條件、表面邊界條件、底部邊界條件、干-濕邊界條件。

3.3.1 開邊界條件

本次模擬模型的開邊界條件為水域邊界條件，即管道路由區(qū)所處的珠江口水域。在此邊界上，給定流速數(shù)據(jù)或給定潮位數(shù)據(jù)。本研究中，小區(qū)域模型的水域邊界條件由經(jīng)過數(shù)據(jù)驗證后的大區(qū)域模型計算得到。

3.3.2 閉邊界條件

本次模擬模型的閉邊界條件為水-陸邊界條件，即海底管道所處的海域與深圳近岸陸地組成的邊界條件。在此邊界上，水質(zhì)點的法向流速為0，即：Vn=0。

3.3.3 表面邊界條件

本次模擬模型的表面邊界條件代表風(fēng)場、氣壓場對水體自由表面的影響，本模型中采用已有的風(fēng)應(yīng)力參數(shù)化方法，對低風(fēng)速、中風(fēng)速及高風(fēng)速條件下的風(fēng)場數(shù)據(jù)進行擬合，表面邊界條件的表達式如下：

式中：Cd為風(fēng)力拖曳系數(shù)；w10為海平面上10 m 處的風(fēng)速，m/s；ca=0.001 255、cb=0.002 425、wa=7 m/s、wb=25 m/s為經(jīng)驗系數(shù)。

3.3.4 底部邊界條件

本次模擬模型中，底部邊界條件通過輸入曼寧系數(shù)M確定，表達式如下：

3.3.5 干-濕邊界條件

本次模擬模型中，對干-濕邊界條件的處理采用的是動邊界方法，在模型計算過程中，模型系統(tǒng)會監(jiān)視每一個網(wǎng)格單元的水深變化值，依據(jù)對干邊界（Dry）、漫水區(qū)（Flood）及濕水區(qū)（Wet）預(yù)先設(shè)定的不同水深值，實時判斷出計算單元的水深類型，并采取相應(yīng)的處理方法[16,17,20-22]。如果監(jiān)測到網(wǎng)格單元的水深值小于干邊界值，則系統(tǒng)將把該網(wǎng)格單元從計算中移除，輸入該單元的動量通量為0[21-25]。

4 工程實測驗證

為了驗證上述海底管道路由區(qū)水動力沖淤數(shù)值模擬模型的有效性，將模型的數(shù)值計算結(jié)果與采用多波束測深儀獲取的工程實測結(jié)果進行對比。

依據(jù)靠近海底管道地理位置附近的赤灣站和珠海站歷年的水文資料統(tǒng)計分析結(jié)果，參考JTS/T 231—2—2010《海岸與河口潮流泥沙數(shù)值模擬技術(shù)規(guī)程》篩掉小波高數(shù)據(jù)，將代表性浪向波高加權(quán)得到研究海域的年均波浪場數(shù)據(jù)[24-26]。

正常海況下（以2019 年對管道路由區(qū)的現(xiàn)場潮位、潮流、懸沙、表層沉積物調(diào)查數(shù)據(jù)以及赤灣站、珠海站的年均波浪、風(fēng)速等數(shù)據(jù)為基礎(chǔ)），采砂區(qū)與管道路由區(qū)在現(xiàn)狀水深條件下1年的水動力沖淤特征模擬模型計算結(jié)果顯示：

采砂區(qū)形成后，對周邊海域沖淤環(huán)境的影響主要體現(xiàn)為采砂坑的南、北兩側(cè)以及東側(cè)為沖刷環(huán)境，海底局部沖刷強度約為1 m/a，采砂坑內(nèi)為淤積環(huán)境，海底最大淤積強度達1 m/a以上。

采砂區(qū)形成后，在水動力作用下，采砂坑北側(cè)邊界向西北方向擴展，遠(yuǎn)離管道走向，對管道安全無影響。采砂區(qū)南側(cè)邊界向東南方向擴展，與管道走向有一定重疊，導(dǎo)致銅鼓水道附近KP 12-KP 14 段管道路由區(qū)海底侵蝕加劇，該段侵蝕強度達到約0.2 m/a，其他管道段仍呈現(xiàn)淤積特征。需特別注意的是，與管道中部平行的、新形成的采砂區(qū)東側(cè)邊界形態(tài)較為曲折，與該區(qū)域的漲、落潮主流向斜交，導(dǎo)致水動力沖淤作用會對該段采砂區(qū)邊界的突出部分進行沖淤調(diào)整，引起凸體局部沖刷，從而使得采砂區(qū)東側(cè)邊界向靠近管道方向擴展。

沖刷1 年后，采砂區(qū)北側(cè)邊界呈環(huán)狀向西北方向擴展1 km，采砂區(qū)南側(cè)邊界呈環(huán)狀向東南方向擴展1 km，采砂區(qū)東側(cè)邊界向靠近管道方向擴展40 m，擴展段為與管道KP 4+120-KP 11+120段之間平行的部分，距管道最近段為KP 5+300-KP 5+600 之間的部分。正常海況下管道路由區(qū)模擬計算1年周期后的水動力沖淤特征如圖6所示。

圖6 基于數(shù)值模擬模型計算的管道路由區(qū)水動力沖淤特征圖（圖中紅線為管道位置）

將2019—2020 年兩期次、間隔1 年的管道路由區(qū)海底地形實測結(jié)果與數(shù)值模擬模型計算1年的結(jié)果進行對比驗證顯示，實測數(shù)據(jù)與模擬值具有很高的符合程度，數(shù)值模擬的預(yù)測準(zhǔn)確度達95%以上，表明數(shù)值模型有效。實測數(shù)據(jù)與數(shù)值模擬值對比見表2。

表2 管道路由區(qū)部分典型位置的沖淤變化實測值與數(shù)值模擬值對比

5 結(jié)論

針對深圳沿海某海底管道路由區(qū)的海床破壞引起海底沖刷、水動力侵蝕變化以及沖淤特征演化等影響管道路由區(qū)安全的問題，研究了適用于該海底管道路由區(qū)特征的水動力沖淤特征數(shù)值模擬方法，通過模型模擬了正常海況下一年周期的海底沖淤變化，并根據(jù)2019—2020 年兩期次、間隔1年的海底地形沖淤變化實測數(shù)據(jù)對數(shù)值模擬模型進行了驗證，得出以下結(jié)論：

1）研究提出了基于水動力數(shù)值模型分析的路由區(qū)沖淤特性計算方法，設(shè)計了二維水流泥沙數(shù)學(xué)模型，實現(xiàn)了符合管道路由區(qū)工程實際的數(shù)值模擬計算。

2）利用管道路由區(qū)實測的水深地形、潮流、懸沙、表層沉積物組成數(shù)據(jù)確定的不同尺度數(shù)值模型的邊界條件在模擬計算中有效，模擬了不同模型尺度下管道路由區(qū)的沖淤狀況，路由區(qū)現(xiàn)場實測海底地形變化數(shù)據(jù)驗證了模型方法的正確性與可靠性。

3）形成的技術(shù)方法可解決海底管道路由區(qū)不同尺度范圍內(nèi)水動力沖淤動態(tài)變化影響的分析難題，可為海底管道保護與運維評估提供技術(shù)支持。