張芝永，劉光生，曾劍

(1.浙江省水利河口研究院浙江省河口海岸重點(diǎn)實(shí)驗(yàn)室，浙江杭州310020;2.浙江大學(xué)建筑工程學(xué)院，浙江杭州310058)

海底管線作為海洋油氣開發(fā)的重要組成部分，是長(zhǎng)距離輸送油氣最經(jīng)濟(jì)的方式之一，它被稱為“海上生命線”。海底管線附近的水動(dòng)力條件引起的海底管線局部沖刷是造成管線事故的一個(gè)主要原因。目前，許多學(xué)者對(duì)于海底管道的局部沖刷問(wèn)題做了大量詳細(xì)、系統(tǒng)的研究，但這些研究所針對(duì)的問(wèn)題主要集中在海底管道局部最大沖刷平衡深度［1-4］及局部沖刷發(fā)展歷程［5-10］方面。在沖刷機(jī)理方面，Chiew［11］應(yīng)用物理模型試驗(yàn)對(duì)發(fā)現(xiàn)管涌是引起管線沖刷的主要因素，Sumer等［12］根據(jù)管涌發(fā)生的臨界條件，通過(guò)物理試驗(yàn)得到了水流、波浪作用下的管線沖刷臨界條件。Zang等［13］利用數(shù)值手段對(duì)波浪、水流作用下管線沖刷的臨界條件進(jìn)行探討。Zhang等［14］基于多孔介質(zhì)假設(shè)，建立了海床滲流場(chǎng)數(shù)值模型，通過(guò)與前人試驗(yàn)結(jié)果對(duì)比驗(yàn)證了該方法的可靠性。

綜上所述，雖然文獻(xiàn)［13］基于數(shù)值模擬方法對(duì)波浪作用下沖刷臨界條件進(jìn)行了研究。但文獻(xiàn)中提到的數(shù)值模擬方法需要分別求解水動(dòng)力場(chǎng)和滲流場(chǎng)，整個(gè)過(guò)程較為繁瑣，而且其未考慮水面變化和垂向流速的影響，與實(shí)際情況有差別。鑒于此，本文提出了一種聯(lián)合求解水動(dòng)力場(chǎng)－滲流場(chǎng)的新方法，該方法將海床假設(shè)為多孔介質(zhì)，然后通過(guò)對(duì)控制方程N(yùn)-S方程中添加源項(xiàng)來(lái)實(shí)現(xiàn)海床沙顆粒對(duì)流體的阻滯作用，這樣水動(dòng)力場(chǎng)和滲流場(chǎng)就可以通過(guò)求解形式大致相同的控制方程來(lái)得到。在此基礎(chǔ)上，根據(jù)線性波理論和VOF方法，對(duì)波浪作用下海底管線周圍水動(dòng)力場(chǎng)和滲流場(chǎng)進(jìn)行耦合求解。

1 數(shù)學(xué)模型

1.1 水動(dòng)力控制方程

無(wú)論是水動(dòng)力場(chǎng)還是滲流場(chǎng)，其控制方程均為雷諾平均N-S方程，包括連續(xù)性方程和動(dòng)量方程。連續(xù)性方程:

動(dòng)量方程:

其中:

式中:u為流體速度;i，j取值1，2;ρ是流體的密度; P是作用在流體微元上的壓力;ν是流體的運(yùn)動(dòng)粘滯系數(shù);是i方向速度脈動(dòng)值;為雷諾應(yīng)力張量;νT紊流運(yùn)動(dòng)粘滯系數(shù);k為紊動(dòng)能;δij是克羅內(nèi)克爾符號(hào)。

1.2 湍流模型

紊流模型采用標(biāo)準(zhǔn)k-ε模型，k方程:

ε方程:

式中:ε為耗散率;Gk是由于平均速度梯度引起的湍動(dòng)能k的產(chǎn)生項(xiàng);μ是流體粘度;μt是紊動(dòng)粘度; C1ε和C2ε為經(jīng)驗(yàn)常數(shù)，取C1ε=1.44，C2ε=1.92;σk和σε分別是與湍動(dòng)能k和耗散率ε對(duì)應(yīng)的Prandtl數(shù)，取σk=1.0，σε=1.3。

1.3 多孔介質(zhì)模型

對(duì)于海床區(qū)域內(nèi)滲流場(chǎng)的求解，可將海床視為均勻、各向同性的多孔介質(zhì)，在動(dòng)量方程(2)添加1個(gè)動(dòng)量阻滯源項(xiàng)來(lái)實(shí)現(xiàn)海床泥沙顆粒的阻滯作用，該動(dòng)量源項(xiàng)公式為

式中:等式右邊括號(hào)內(nèi)第1項(xiàng)為粘性損失，α為粘滯阻力系數(shù);等式右邊括號(hào)內(nèi)第2項(xiàng)為慣性損失。由于海床中滲流流速都是比較小的，其慣性損失也較小，因此在這里只考慮粘滯損失，忽略慣性損失。粘滯阻力系數(shù)α定義為

式中:dp為泥沙平均粒徑，在這里用中值粒徑d50代替;n為孔隙率。

1.4 邊界造波法

對(duì)于波浪模擬，本文通過(guò)給定入口邊界的流速及波高時(shí)程變化來(lái)進(jìn)行造波，根據(jù)線性波理論，入射邊界處的水平速度、垂向速度和水深可以表示為

式中:ω=2π/T，A、k、h、T分別為波浪波高、波數(shù)、靜水深和周期。

1.5 消波邊界

對(duì)于數(shù)值模擬來(lái)說(shuō)，不可能把整個(gè)海域都作為計(jì)算區(qū)域，需要在計(jì)算域末端將波能消去，以避免反射回的波浪影響計(jì)算結(jié)果的準(zhǔn)確性。在數(shù)值波浪水槽中可采用阻尼層吸收邊界條件處理開邊界，使波浪在阻尼層內(nèi)衰減，從而消除反射波。在這里阻尼層靠近下游出口位置，阻尼層的長(zhǎng)度約為2～3倍的波長(zhǎng)。阻尼層內(nèi)動(dòng)量方程中源項(xiàng)和衰減系數(shù)分別取如下形式:

式中:x0為阻尼層的起點(diǎn)橫坐標(biāo)，x1為阻尼層的終點(diǎn)橫坐標(biāo)，xL為阻尼層內(nèi)任意點(diǎn)的橫坐標(biāo)。

2 水動(dòng)力場(chǎng)－滲流場(chǎng)耦合模型驗(yàn)證

2.1 模型驗(yàn)證1

本文首先對(duì)水流作用下管線周圍水動(dòng)力場(chǎng)－滲流場(chǎng)進(jìn)行聯(lián)合數(shù)值求解來(lái)驗(yàn)證本耦合模型的可靠性，然后在此基礎(chǔ)上將水流條件改為波浪條件，對(duì)波浪作用下的水動(dòng)力場(chǎng)－滲流場(chǎng)進(jìn)行驗(yàn)證。

為便于對(duì)比，計(jì)算中采用無(wú)量綱的壓力系數(shù)，其公式為

式中:p為計(jì)算區(qū)域內(nèi)任意一點(diǎn)壓力;p0為參考?jí)毫?，在這里取初始靜水壓力;u0為速度入口水面位置處流速，波浪情況時(shí)將u0改為近底最大流速umax。

Chiew應(yīng)用物理模型試驗(yàn)對(duì)海底管線沖刷機(jī)理進(jìn)行了研究并對(duì)部分埋入海床的海底管線表面壓分布進(jìn)行了測(cè)量。在這里選取其中一個(gè)工況來(lái)對(duì)本文的水動(dòng)力場(chǎng)－滲流場(chǎng)耦合模型進(jìn)行驗(yàn)證。該工況為來(lái)流流速為 0.4 m/s，水深為 0.3 m，管線直徑0.05 m，相對(duì)水深為6.0，水面變化不大，因此在這里不考慮水面的變化。其管線埋深為0.5D，試驗(yàn)?zāi)嗌沉綖?.7 mm，密度ρs取2 600 kg/m3，孔隙率n= 0.43，根據(jù)太沙基公式，發(fā)生管涌流土的臨界水力梯度為

式中:s為泥沙比重，s=ρs/ρ。將參數(shù)代入上式得到海床的臨界水力梯度ic約為0.91。

上游入口邊界和下游邊界到管線中心位置的距離均為20D，水流計(jì)算區(qū)域水深取6D，滲流計(jì)算區(qū)域豎直方向距離取8D，整個(gè)計(jì)算區(qū)域采用結(jié)構(gòu)化網(wǎng)格進(jìn)行離散。

控制方程采用有限體積法進(jìn)行離散，壓力和速度采用SIMPLE算法進(jìn)行耦合求解，在水流區(qū)域內(nèi)求解標(biāo)準(zhǔn)N-S控制方程的同時(shí)，海床區(qū)域內(nèi)滲流場(chǎng)通過(guò)求解添加了阻滯源項(xiàng)的控制方程。

圖1為管線周圍海床內(nèi)外流場(chǎng)圖，以管線最上游端為0°起始位置，順時(shí)針方向逐漸增大。這樣管線最頂端位置為90°，管線最下游端位置為180°。從圖中可以看出，在海床以上，水流自由流動(dòng)，由于管線的阻擋發(fā)生繞流現(xiàn)象，在管線后方出現(xiàn)一較大的渦，而在管線前方根部同樣出現(xiàn)一個(gè)較小渦體。而在海床內(nèi)，由于上下游壓差的存在，滲流由管道上游流向管道下游。而且由于阻滯作用速度較小，并未出現(xiàn)渦旋。圖2為管線表面壓力分布情況，可以看出計(jì)算結(jié)果與試驗(yàn)結(jié)果的規(guī)律較為一致，在裸露管線部分(0°～180°)內(nèi)，其最大壓力位于海床以上管線迎流面范圍(0°～45°)，在接近90°位置，即管線最頂端位置，壓力最小。在掩埋管線部分(180°～360°)內(nèi)，管線表面壓力分布與裸露部分壓力分布完全不同，表面壓力沿著360°～180°的弧線方向逐漸減小。而且從曲線上可以看出兩端點(diǎn)位置即180°和360°位置的壓力變化率要大于掩埋管線部分的其他的位置。這說(shuō)明兩端點(diǎn)處的水力梯度要大于其他位置的水力梯度，通過(guò)求解，下游端滲流出口處水力梯度為0.27，明顯小于管涌發(fā)生的臨界水力梯度0.91，因而在掩埋0.5D情況下，沖刷不會(huì)發(fā)生。這一結(jié)論在試驗(yàn)中得到了驗(yàn)證。

圖1 管線周圍流場(chǎng)圖Fig.1 Flow vectors around pipeline

圖2 管道表面壓力分布情況Fig.2 Pressure distributions along the pipeline surface

2.2 模型驗(yàn)證2

應(yīng)用Sumer的試驗(yàn)來(lái)進(jìn)行波浪作用下的水動(dòng)力場(chǎng)－滲流場(chǎng)模擬驗(yàn)證，具體參數(shù):管線直徑D為0.1 m，管線埋深0.064D，波高A為0.17 m，周期T為4 s，水深0.33 m，整個(gè)計(jì)算區(qū)域長(zhǎng)度為6倍的波長(zhǎng)，管線位于中間位置。計(jì)算區(qū)域采用結(jié)構(gòu)化網(wǎng)格進(jìn)行離散。

圖3為管線上下游端的壓差變化情況，可以看出數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果基本吻合，在波浪作用下管線兩端壓差周期性變化，在正壓差情況下，滲流由管道上游向下游流動(dòng)，負(fù)壓差情況下，滲流則由管道下游向上游流動(dòng)，由圖3發(fā)現(xiàn)，正壓差幅值要大于負(fù)壓差幅值，因此取壓差為最大幅值時(shí)的管線表面壓力分布，然后依據(jù)水流條件下的水力梯度求解過(guò)程即可得到該波浪工況下的最大水力梯度值。當(dāng)這一水力梯度小于管涌發(fā)生的臨界水力梯度時(shí)，沖刷不會(huì)發(fā)生。

圖3 上下游端壓差變化情況Fig.3 Pressure difference between upstream and downstream point

3 結(jié)果分析

應(yīng)用上述模型對(duì)不同波浪、不同埋深情況下的管線周圍水動(dòng)力場(chǎng)－滲流場(chǎng)進(jìn)行耦合求解。

3.1 波浪KC數(shù)對(duì)壓差影響

對(duì)于波浪來(lái)說(shuō)，常常用無(wú)量綱KC(Keulegan Carpenter)數(shù)來(lái)表示其特征要素，其定義為

式中:umax為近底最大流速，其雷諾數(shù)Re=umaxD/ν，范圍在16 699～45 058，屬于紊流范圍。

圖4 壓差系數(shù)隨KC變化情況Fig.4 Variations of pressure difference coefficient with KC

3.2 相對(duì)埋深e/D對(duì)壓差影響

埋深不同的管線兩端壓差系數(shù)大小也有所不同。圖5為不同埋深情況下壓差系數(shù)的變化規(guī)律?？梢钥闯?，相同KC下，隨著管道埋深的加大，裸露在水流中的管道高度越來(lái)越小，其對(duì)水體的阻流作用也越來(lái)越小，因此其壓差系數(shù)隨著埋深的加大而逐漸減小。而且呈現(xiàn)近似線性遞減關(guān)系。

圖5 壓差系數(shù)隨相對(duì)埋深e/D變化情況Fig.5 Variations of pressure difference coefficient with relative buried depth e/D

3.3 水力梯度分析

當(dāng)管線兩端壓差達(dá)到最大時(shí)，依照恒定水流條件下滲流出口處水力梯度求解方法，對(duì)管線表面壓力沿著弧長(zhǎng)進(jìn)行求導(dǎo)，然后通過(guò)換算即可得到各工況下滲流出口處最大水力梯度，圖6列出了KC=11時(shí)，不同埋深情況下滲流出口處最大水力梯度的變化情況。從圖6中可以看出，隨著埋深的加大，由于兩端壓差的減小和滲流路徑的加長(zhǎng)，其水力梯度急劇減小。圖6中虛線為管涌滲透破壞臨界水力梯度曲線，該曲線與水力梯度－埋深關(guān)系曲線交于一點(diǎn)，在該點(diǎn)處，其最大水力梯度正好等于滲透破壞臨界水力梯度，因而此點(diǎn)即為管涌乃至沖刷的臨界點(diǎn)，在該點(diǎn)左側(cè)，其水力梯度大于臨界水力梯度，因而沖刷均發(fā)生，而在右側(cè)則小于臨界水力梯度，沖刷并未發(fā)生。這樣通過(guò)畫出不同KC情況下的最大水力梯度隨埋深變化曲線，并結(jié)合管涌臨界水力梯度值，即可通過(guò)這種圖解方法得到其臨界條件。

圖6 滲流出口處水力梯度隨埋深e/D變化情況Fig.6 Variations of hydraulic gradient with e/D

3.4 沖刷臨界條件

利用3.3節(jié)中方法對(duì)不同KC時(shí)的管線沖刷的臨界條件進(jìn)行了研究。臨界條件不僅與埋深e/D、KC有關(guān)還與泥沙性質(zhì)有關(guān)。Sumer等定義了與泥沙性質(zhì)相關(guān)的一個(gè)無(wú)量綱臨界流速:

在這里同樣采用該無(wú)量綱臨界流速來(lái)進(jìn)行分析，圖7中列出了不同情況下無(wú)量綱臨界流速隨KC和埋深的變化情況?？梢钥闯鰯?shù)值計(jì)算結(jié)果與Sumer等［12］的試驗(yàn)結(jié)果基本吻合。由前述分析可知，管線埋深越深，管線對(duì)流體的阻礙作用越弱，管線兩端壓差越小，進(jìn)而導(dǎo)致滲流出口處水力梯度較小，在這種情況下要引起沖刷就必然增大水動(dòng)力條件。因此，隨著管線埋深的加大，其無(wú)量綱臨界流速也在逐漸增大。對(duì)于固定埋深的情況，隨著KC的增大，壓差逐漸減小，因而其臨界無(wú)量綱流速也是需要逐漸增大，Sumer指出當(dāng)KC趨向很大時(shí)，其沖刷臨界流速趨于恒定流條件下的沖刷臨界流速。對(duì)于恒定流條件下的臨界流速，其經(jīng)驗(yàn)公式為

而圖7中不同KC情況下的臨界流速隨KC近似線性變化，因此在這里利用KC對(duì)式(17)進(jìn)行擬合修正得到波浪作用下無(wú)量綱臨界流速的近似計(jì)算公式:

通過(guò)圖7可以看出擬合公式的計(jì)算結(jié)果與數(shù)值計(jì)算結(jié)果比較接近，而且反映出了無(wú)量綱臨界流速隨KC和埋深的變化規(guī)律。這說(shuō)明可以應(yīng)用此公式來(lái)進(jìn)行波浪作用下沖刷臨界條件的計(jì)算。

4 結(jié)論

基于多孔介質(zhì)假設(shè)和邊界造波法，建立了波浪作用下海底管線周圍水動(dòng)力場(chǎng)－滲流場(chǎng)耦合數(shù)值模型。通過(guò)與前人試驗(yàn)結(jié)果對(duì)比驗(yàn)證了耦合模型的可靠性。在此基礎(chǔ)上，應(yīng)用該模型波浪作用下海底管線沖刷臨界條件進(jìn)行了研究，結(jié)果表明:

1)當(dāng)KC＜20時(shí)，隨著KC的增大，管線兩端壓差系數(shù)急劇減小，而在KC≥20后，壓差隨KC的變化幅度越來(lái)越小并逐漸趨于穩(wěn)定;

2)相同KC情況下，隨著埋深的加大，壓差系數(shù)近似線性遞減，而滲流出口處水力梯度則明顯減小。

3)波浪作用下管線沖刷無(wú)量綱臨界流速隨著KC和埋深的增加而均增大，其值可通過(guò)式(18)近似求得，這為實(shí)際工程中海底管線防護(hù)決策提供了技術(shù)支持。

［1］SUMER B M，JENSEN H R，F(xiàn)REDSE J.Effect of lee-wake on scour below pipelines in current［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，1988，114 (5):599-614.

［2］MONCADA-M A T，AGUIRRE-PE J.Scour below pipeline in river crossings［J］.Journal of Hydraulic Engineering，1999，125(9):953-958.

［3］DEY S，SINGH N P.Clear-water scour depth below underwater pipelines［J］.Journal of Hydro-environment Research，2007，1(2):157-162.

［4］李寧，谷凡.海底管線局部沖刷平衡深度Chao解析法的改進(jìn)［J］.大連海事大學(xué)學(xué)報(bào)，2012，38(1):77-80.

LI Ning，GU Fan.Improvement of Chao-Hennesy analytical method for estimating the equilibrium depth of local scour around submarine pipeline［J］.Journal of Dalian Maritime University，2012，38(1):77-80.

［5］BRRS B.Numerical modeling of flow and scour at pipelines［J］.Journal of Hydraulic Engineering，1999，125(5): 511-523.

［6］LI Fangjun，CHENG Liang.Numerical model for local scour under offshore pipelines［J］.Journal of Hydraulic Engineering，1999，125(4):400-406.

［7］LIANG Dongfang，CHENG Liang，LI Fangjun.Numerical modeling of scour below a pipeline in currents:partⅡ.Scour simulation［J］.Coastal Engineering，2004，52(1): 43-62.

［8］GAO Fuping，YANG Bing，WU Yingxiang，et al.Steady current induced seabed scour around a vibrating pipeline［J］.Applied Ocean Research，2006，28(5):291-298.

［9］文君鋒，勾瑩，宋偉華，等.海流作用下海底管道局部沖刷數(shù)值分析［J］.海洋工程，2012，30(1):75-82.

WEN Junfeng，GOU Ying，SONG Weihua，et al.Numerical examination of steady current-induced local scour around submarine pipeline［J］.The Ocean Engineering，2012，30 (1):75-82.

［10］張芝永，拾兵.泥沙局部沖淤二維數(shù)值模擬仿真［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2013，34(2):145-150.

ZHANG Zhiyong，SHI Bing.Two-dimensional numerical simulation on local mud＆sand erosion and deposition［J］.Journal of Harbin Engineering University，2013，34(2): 145-150.

［11］CHIEW Y M.Mechanics of local scour around submarine pipelines［J］.Journal of Hydraulic Engineering，1990，116 (4):515-529.

［12］SUMER B M，TRUELSEN C，SICHMANN T，et al.Onset of scour below pipelines and self-burial［J］.Coastal Engineering，2001，42(4):313-335.

［13］ZANG Zhiping，CHENG Liang，ZHAO Ming，et al.A numerical model for onset of scour below offshore pipelines［J］.Coastal Engineering，2009，56(4):458-466.

［14］ZHANG Zhiyong，SHI Bing，GUO Yakun，et al.Numerical investigation on critical length of impermeable plate below underwater pipeline under steady current［J］.Science China(Technological Sciences)，2013，56(5):1232-1240.

［15］劉莎莎，顧煜炯，惠萬(wàn)馨，等.基于邊界造波法的波浪數(shù)值模擬［J］.可再生能源，2013，31(2):100-103.

LIU Shasha，GU Yujiong，HUI Wanxin，et al.Wave numerical simulation based on wave-generation method of defining inlet boundary conditions［J］.Renewable Energy Resources，2013，31(2):100-103.