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

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

2.中國石油集團(tuán)海洋工程重點(diǎn)實(shí)驗(yàn)室，天津 300451

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

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

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

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

1 研究區(qū)域概況

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

圖1 海底管道地理位置

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

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

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

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

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

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

2.2 數(shù)值模型研究

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

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

式中：t為時(shí)間，s；x、y、z是笛卡爾坐標(biāo)系坐標(biāo)，m；η為水面標(biāo)高，m；p為水體壓強(qiáng)，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。

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

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

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

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

2.3 泥沙輸運(yùn)模型研究

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

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

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

式中：c為泥沙的質(zhì)量濃度，g/m3；ωs為泥沙顆粒在水中的沉降速度，m/s；σc為紊流Schmidt數(shù)，表示泥沙擴(kuò)散系數(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時(shí)，利用式（7）計(jì)算海底邊界的泥沙質(zhì)量濃度；當(dāng)T≤0時(shí)，海底泥沙的質(zhì)量濃度在邊界法線方向梯度為0。在運(yùn)動(dòng)水體的圓柱表面和自由面，水體邊界法線方向的泥沙流量為0，泥沙質(zhì)量濃度計(jì)算方程[7-9,20]如下：

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

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

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

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

本研究區(qū)推移質(zhì)的輸運(yùn)受泥沙密度、流體密度及流體速度控制[20-22]，本次數(shù)值模擬中推移質(zhì)輸沙率[23-26]的計(jì)算公式[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為總泥沙輸運(yùn)率，kg/s，等于推移質(zhì)輸沙率qb和懸移質(zhì)輸沙率qs之和[18-22]。

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

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

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

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

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

圖4 大區(qū)域模型計(jì)算范圍與網(wǎng)格劃分圖（黑色區(qū)域?yàn)檠芯繀^(qū)）

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

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

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

3.3 邊界條件設(shè)置

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

3.3.1 開邊界條件

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

3.3.2 閉邊界條件

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

3.3.3 表面邊界條件

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

式中：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)驗(yàn)系數(shù)。

3.3.4 底部邊界條件

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

3.3.5 干-濕邊界條件

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

4 工程實(shí)測驗(yàn)證

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

依據(jù)靠近海底管道地理位置附近的赤灣站和珠海站歷年的水文資料統(tǒng)計(jì)分析結(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年的水動(dòng)力沖淤特征模擬模型計(jì)算結(jié)果顯示：

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

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

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

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

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

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

5 結(jié)論

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

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

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

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