夏華永,王俊勤

(1.國家海洋局南海預(yù)報(bào)中心,廣東廣州510300; 2.中海石油(中國)有限公司北京研究中心,北京100027)

波浪作用下管道沖刷暴發(fā)計(jì)算的邊界元方法

夏華永1,王俊勤2

(1.國家海洋局南海預(yù)報(bào)中心,廣東廣州510300; 2.中海石油(中國)有限公司北京研究中心,北京100027)

管道沖刷暴發(fā)臨界條件是海底管道設(shè)計(jì)和運(yùn)營的重要參數(shù)。Sumer等給出了3個(gè)Keulegan-Carpenter (KC)數(shù)的波浪作用下沖刷暴發(fā)的試驗(yàn)曲線。但對于任意KC數(shù),Sumer等的方法并不準(zhǔn)確。因此,建立波浪作用下沖刷暴發(fā)條件計(jì)算方法,具有重要意義。因管道與床面相交形成了復(fù)雜的幾何形狀,管道與床面的交叉點(diǎn)是數(shù)值奇點(diǎn),因而管道沖刷暴發(fā)臨界條件計(jì)算較為復(fù)雜。本研究采用邊界元方法計(jì)算波浪運(yùn)動(dòng)及滲流壓力。邊界元可準(zhǔn)確地?cái)M合海底與管壁邊界,方便地處理管道與海床交叉點(diǎn)的數(shù)值奇點(diǎn)問題,并可直接計(jì)算出奇點(diǎn)處的垂向壓力梯度,準(zhǔn)確地計(jì)算沖刷暴發(fā)的臨界條件。本研究計(jì)算結(jié)果與實(shí)測值符合良好。

海底管道;沖刷暴發(fā);邊界元;數(shù)值奇點(diǎn)

海底管道懸空是管道運(yùn)行極大的安全隱患,因此,管道局部沖刷是管道設(shè)計(jì)與運(yùn)行必需掌握的參數(shù)。目前,已發(fā)展了多種海底管道沖刷模型,如勢流模型、大渦模擬模型、紊流模型、渦度模型、兩相流模型及格點(diǎn)波爾茲曼方法模型[1-6]。模擬過程中,這些模型都假定管道下面與海床之間已有一個(gè)小的縫隙,即假定了管道沖刷已經(jīng)發(fā)生,然后模擬在水流與波浪作用下的海床沖刷過程及平衡剖面形態(tài)。但是,鋪設(shè)在海底的管道在其自重作用下會發(fā)生一定程度的自埋,管道是否會發(fā)生沖刷呢?這需要首先計(jì)算其沖刷暴發(fā)的臨界條件。

目前,對于管道沖刷暴發(fā)的機(jī)制與臨界條件,已有多種試驗(yàn)研究[7-11]。試驗(yàn)結(jié)果已揭示了沖刷暴發(fā)的機(jī)制,確定了沖刷暴發(fā)的判別條件,建立了沖刷暴發(fā)的計(jì)算方法。試驗(yàn)表明,管道上游渦旋與管道下面滲流的聯(lián)合作用是沖刷暴發(fā)的原因,后者起主要作用。管道與海床的下游交叉點(diǎn)處,滲流垂向壓力梯度大于泥沙浮重時(shí),產(chǎn)生管涌,沖刷暴發(fā)。Sumer等[11]根據(jù)試驗(yàn)資料,整理了沖刷暴發(fā)臨界條件的計(jì)算方法,單向流作用下,給出了沖刷暴發(fā)臨界條件的經(jīng)驗(yàn)公式,計(jì)算方便,對于波浪作用下沖刷暴發(fā)的臨界條件,則給出了3個(gè)Keulegan-Carpenter(KC)數(shù)的圖示曲線,可根據(jù)圖示曲線得到臨界條件。

對于沖刷暴發(fā)臨界條件的模型計(jì)算,就本文作者所知,只有極少的研究成果。楊兵等[12]模擬分析了管道周圍水體及滲流的壓力分布,重現(xiàn)了管道沖刷的暴發(fā)機(jī)制。Liang和Cheng[13],Zang等[14]建立模型,模擬計(jì)算了單向流與波浪作用下的管道沖刷暴發(fā)的臨界條件。管道輕微埋在海床上,在管道與海床的交叉點(diǎn)形成了復(fù)雜的幾何形狀。沖刷暴發(fā)計(jì)算需要處理好這一復(fù)雜的邊界,而且處理好這些邊界對流場及壓力模擬非常重要,特別是在交叉點(diǎn)。結(jié)構(gòu)網(wǎng)格很難處理好這一幾何形狀。非結(jié)構(gòu)網(wǎng)格可以處理這一邊界。但是,在滲流計(jì)算中,交叉點(diǎn)是滲流壓力梯度的奇異點(diǎn),差分法與有限元法都難于計(jì)算交叉點(diǎn)的垂向壓力梯度,因而也影響沖刷暴發(fā)臨界的預(yù)測精度。Liang和Cheng[13],Zang等[14]都通過建立下游交叉點(diǎn)處滲流垂向壓力梯度與兩個(gè)交叉點(diǎn)之間平均壓力差之間的關(guān)系,進(jìn)一步計(jì)算沖刷暴發(fā)的臨界條件。

目前,沖刷計(jì)算一般采用Navier-Stokes方程,嵌入紊流封閉子模型提供參數(shù)化方案。滲流的控制方程則為Laplace方程。但是,波浪運(yùn)動(dòng)通常當(dāng)作勢流來處理。勢流與滲流都可以準(zhǔn)確地用邊界元方法計(jì)算。對于Laplace方程,邊界元法是一種非常有效的計(jì)算方法[15]。邊界元具有準(zhǔn)確擬合邊界,降低計(jì)算維數(shù),大幅度的減小計(jì)算量、可得到精確解、達(dá)到任意分辨率等優(yōu)勢。邊界元法可以處理奇點(diǎn)問題,給出交叉點(diǎn)準(zhǔn)確的滲流垂向壓力梯度。對于任意KC數(shù)下臨界沖刷條件,Sumer等[11]建立的圖示并不能方便的給出。建立簡單而準(zhǔn)確的波浪作用下沖刷暴發(fā)臨界條件的計(jì)算方法,有其明確的應(yīng)用價(jià)值。

(李 燕 編輯)

1 沖刷暴發(fā)機(jī)制

Sumer等[11]對波流作用下的管道沖刷暴發(fā)機(jī)制已做了全面的總結(jié),為了論文完整性,對沖刷暴發(fā)臨界條件在此簡單介紹一下。海底鋪設(shè)了管道后,改變了局部的流態(tài),流態(tài)的改變導(dǎo)致管道附近水體及床面上的壓力變化。床面的壓力分布不均勻及管道對滲流結(jié)構(gòu)的影響導(dǎo)致了管道下沙體較大的壓力梯度存在(圖1)。圖1中,管道直徑為D,管道埋入海床的深度為e,B點(diǎn)與A點(diǎn)分別為管道與海床的上、下游交叉點(diǎn)。當(dāng)在管道下游管壁與床面的交叉點(diǎn)(圖1中的A點(diǎn))的垂向壓力梯度大于沙土水下重量時(shí),管道沖刷暴發(fā)

式中,PA為交叉點(diǎn)A處滲流壓強(qiáng);s0=ρs/ρ,ρ為海水密度,ρs泥沙密度;n為沙體的孔隙率;g為重力加速度;W為泥沙的水下重量。

Zang等[14]取滲流出流點(diǎn)A處的垂向壓力梯度與管道下方的平均壓強(qiáng)相關(guān)

圖1 管道下滲流分布示意圖Fig.1 A sketch of seepage flow distribution beneath the marine pipeline

式中,λA是一個(gè)校正系數(shù);s是沿管壁距離。試驗(yàn)結(jié)果顯示,管道下滲流壓強(qiáng)變化劇烈,因此,λA是一個(gè)變化復(fù)雜的系數(shù)。取管壁下平均壓強(qiáng)為

式(3)中,D為管道直徑;u0為海底未受擾動(dòng)時(shí)的流速;ΔCP0為A,B兩點(diǎn)之間的壓力降幅系數(shù),CPA,CPB為壓力系數(shù)。由式(1)～式(3),得到?jīng)_刷暴發(fā)條件的替代表述

在波浪作用下,管道沖刷暴發(fā)與Keulegan-Carpenter數(shù)有關(guān),也與管道的埋設(shè)深度有關(guān),沖刷暴發(fā)條件為

式(4)中,um為床面波浪質(zhì)點(diǎn)軌跡的最大運(yùn)動(dòng)速度;KC為Keulegan-Carpenter數(shù),KC=umT/D,T為波浪周期。

2 控制方程

假定波浪傳播方向垂直于管道路由。波浪運(yùn)動(dòng)是無旋的,存在流速勢函數(shù),其控制方程為Laplace方程

對于微幅波,波浪勢為

假定滲流符合Darcy定律,滲流控制方程為

式中,P(x,y)為壓強(qiáng)。在平坦海床上,波浪動(dòng)水壓強(qiáng)為

顯然,式(8)是式(7)的一個(gè)邊界。根據(jù)邊界條件的形式,將方程(7)的解表示為P=f(y)cos(kx-ωt),代入到式(7)中,并進(jìn)一步處理,得到形式如下的解P=c1exp( ky)cos( kx-ωt)+c2exp(-ky)cos(kx-ωt),當(dāng)沙層無限厚時(shí),解為

如果在y=-d,存在不滲水層,即?P/?y=0,則解為

3 邊界元法

本文計(jì)算區(qū)域如圖2所示。波浪計(jì)算的邊界為Γ1,Γ2,Γ3,Γ4。Γ1由海底與未埋藏的管壁構(gòu)成。在波浪計(jì)算中,假定管道不影響較遠(yuǎn)水體中的質(zhì)點(diǎn)運(yùn)動(dòng)規(guī)律,邊界條件取為

在上式中,n表示邊界的外法線方向。上述邊界條件全為Neumann型,對于Laplace方程而言,是沒有唯一解的,因此,在波浪的入射斷面上,改為根據(jù)式(6)給定勢函數(shù)值。

圖2 波浪、滲流計(jì)算區(qū)域示意圖Fig.2 A sketch of the computation domain for the sea wave and the seepage flow

圖3 邊界積分單元定義Fig.3 Definition of a boundary integral element

24定?P/?n,在邊界L1上,如果沙層足夠厚,則根據(jù)式(9)給定壓強(qiáng),如果存在不滲水層,則取?P/?n=0。對于Laplace方程,其基礎(chǔ)解為

式(12)中,(x,y)與(η,ζ)為計(jì)算域中的點(diǎn),r為點(diǎn)(x,y)到點(diǎn)(η,ζ)的距離。通過格林積分公式

得到Laplace方程的解

式(13)中,Ω為計(jì)算域;Γ=Γ1+Γ2+Γ3+Γ4,為計(jì)算域的邊界;公式左邊為計(jì)算域中的面積分,右邊為邊界上的線積分。式(14)中,C為積分系數(shù),C=θ/2π,θ為點(diǎn)(η,ζ)鄰域的積分幅度,如(η,ζ)∈ΩΓ,θ= 2π,(η,ζ)∈Γ且光滑,θ=π。

根據(jù)式(16),如已知邊界上的qi(或φi),就可以確定邊界上φi(或qi),邊界上φi與qi都已知時(shí),則采用(14)計(jì)算區(qū)域內(nèi)任意點(diǎn)的勢函數(shù)。

對于邊界元方法而言,計(jì)算區(qū)域的邊界能被準(zhǔn)確擬合,不存在差分法中對奇點(diǎn)難于處理的問題。在管壁與海床接觸的位置(圖1中的A,B點(diǎn))處,以A點(diǎn)為例,與其相接的2個(gè)邊界元記為li-1與li(圖4),那么在A點(diǎn)處,有2個(gè)外法線方向的壓強(qiáng)梯度,其中垂直于管壁的壓強(qiáng)梯度?Pi,2/?n=0,而A點(diǎn)處的壓強(qiáng)Pi也已經(jīng)通過波浪計(jì)算得到,因此,只須要計(jì)算垂直于床面的壓強(qiáng)梯度?Pi,1/?n,而?Pi,1/?n則用于判別管道沖刷暴發(fā)的臨界條件。

數(shù)值試驗(yàn)的結(jié)果表明,在管壁與交叉點(diǎn)鄰近床面,流速與壓強(qiáng)變化劇烈,需要較高的空間分辨率,本文計(jì)算中,在管壁及與交叉點(diǎn)相距0.5D的床面上,邊界元長度取l=πD/72。

計(jì)算過程中,將一個(gè)波浪周期分為N個(gè)時(shí)間步長,Δt=T/N,在t=i(-1)Δt,i=1,2,…,N+1時(shí)計(jì)算波浪,t=i(-1/2)Δt, i=1,2,…,N時(shí)計(jì)算滲流。設(shè)Γ4為波浪入射邊界。計(jì)算流程如下:

1)采用邊界Γ4以式(6)給定波浪勢函數(shù),在Γ1,Γ2,Γ3按式(11)給定波浪勢函數(shù)法向梯度(流速);

2)用邊界元法計(jì)算t=iΔt,(i+1)Δt時(shí)刻的波浪勢函數(shù)方程(5),得到海底(邊界Γ1)的勢函數(shù)φi,φi+1;

圖4 管道與海床交叉處邊界元的處理Fig.4 Treatment of the boundary elements at the cross points of the pipeline and the seabed

3)取Pi+1/2w=-ρ(φi+1-φi)/Δt,計(jì)算t=(i+1/2)Δt時(shí)刻海底的波浪壓強(qiáng),并作為滲流的邊界條件,在邊界L2,L4上,根據(jù)式(9)或(10)給定?P/?n,在邊界L1上,如透水,則據(jù)式(9)給定壓強(qiáng),如果不滲水,則取?P/?n=0;

4)用邊界元法計(jì)算t=(i+1/2)Δt時(shí)刻的滲流控制方程(7),得到海底(邊界L3)的法向壓力梯度?Pi+1/2/?n。

上述計(jì)算過程中,管道與海底交叉點(diǎn)處的法向壓力梯度?P/?n大于泥沙水下重量時(shí),管道沖刷暴發(fā)。

圖5 Test O16中交叉點(diǎn)壓強(qiáng)梯度過程線[11]Fig.5 Time series of pressure gradient at the cross point in Test O16[11]

4 模型驗(yàn)證

本文采用Sumer等[11]的試驗(yàn)數(shù)據(jù)進(jìn)行模型檢驗(yàn)。其試驗(yàn)水槽長26.5 m,寬0.6 m,試驗(yàn)水深0.33 m,水槽中部有一段0.10 m厚,3 m長的沙層,泥沙粒徑0.18 mm,空隙率0.53。試驗(yàn)條件下,沖刷暴發(fā)的臨界條件為-?(P/γ)/?y≥0.77(γ=ρg,g為重力加速度)。Sumer等的試驗(yàn)中,逐漸增大波高,直到?jīng)_刷暴發(fā)為止,Test O16給出了交叉點(diǎn)的壓強(qiáng)梯度過程(圖5)。Test O16的管徑為10 cm,埋深0.064倍管徑,波浪周期為4 s,波高0.17 m。Test O16條件下,本文計(jì)算的交叉點(diǎn)A處的垂向壓強(qiáng)梯度如圖6所示。計(jì)算的最大水頭梯度為1.29,與試驗(yàn)的最大值相當(dāng)。但是,模擬的壓強(qiáng)梯度變化過程是對稱的,而試驗(yàn)結(jié)果是不對稱的,計(jì)算的最小梯度與試驗(yàn)結(jié)果相差較大。其原因在于試驗(yàn)的水深較小,波高較大,管道影響了的波浪的形狀,造成了波浪的不對稱性,因此,試驗(yàn)中的壓強(qiáng)梯度變化是不對稱的,而模型不能重現(xiàn)試驗(yàn)中的不對稱過程。對于水深較大的情形,管道不影響波浪的對稱性,模擬的結(jié)果將更真實(shí)。Test O17條件下,模擬的速度與壓強(qiáng)水頭分布如圖7,交叉點(diǎn)壓強(qiáng)梯度過程見圖6。模擬的最大壓強(qiáng)梯度為0.83,略大于起動(dòng)臨界條件。在波峰時(shí)刻,管道附近水體水平流動(dòng),水體與沙體中的壓強(qiáng)大致對稱分布,交叉點(diǎn)的壓強(qiáng)梯度最小。在T/ 4時(shí)刻,管道上方水體向下流,埋藏管道下方壓強(qiáng)變化劇烈,交叉點(diǎn)的壓強(qiáng)梯度最大。在波谷與3T/4時(shí)刻,速度與壓強(qiáng)分布與波峰與T/4時(shí)刻正好相反。試驗(yàn)結(jié)果中,交叉點(diǎn)垂向壓強(qiáng)梯度達(dá)到,甚至大于臨界值時(shí),沖刷才暴發(fā)。Sumer等[11]對此給出了解釋,波浪過程中,交叉點(diǎn)泥沙暴露在臨界壓強(qiáng)梯度力的時(shí)間較短,而沖刷暴發(fā)需要一個(gè)過程。因此,不同于單向流作用下的沖刷暴發(fā),波浪作用下,垂向壓強(qiáng)梯度需要大于臨界值。模擬結(jié)果,交叉點(diǎn)的壓強(qiáng)梯度都略大于或顯著大于(如Test O16)臨界值。

圖6 模擬交叉點(diǎn)A處壓強(qiáng)梯度過程Fig.6 Simulation of the time series of pressure gradient at cross point A

圖7 波浪質(zhì)點(diǎn)速度與壓強(qiáng)水頭分布Fig.7 Distributions of the wave particle velocity and the pressure head

5 結(jié) 語

邊界元方法能準(zhǔn)確地?cái)M合海底與管壁邊界,可以處理管道與海床交叉點(diǎn)數(shù)值奇點(diǎn)問題,準(zhǔn)確地計(jì)算波浪與滲流。本文采用邊界元方法計(jì)算了波浪作用下管道沖刷暴發(fā)臨界條件,計(jì)算結(jié)果與試驗(yàn)結(jié)果符合良好。模型方法程序設(shè)計(jì)簡單,計(jì)算精度高,可根據(jù)波浪與管道埋深條件,快速計(jì)算交叉點(diǎn)的垂向壓力梯度,判斷管道是否會發(fā)生沖刷,彌補(bǔ)了根據(jù)試驗(yàn)圖示判別的不足。

當(dāng)然,本文方法也有適用范圍。1)本文方法是采用線性波理論計(jì)算水體中的壓力場及線性波作用下的滲流場,而在波浪破碎帶內(nèi),波面發(fā)生了變形,已不能用線性波來描述,本文方法就已經(jīng)失去了其賴以建立的條件,可能產(chǎn)生較大的誤差。對于近岸波浪,非線性波理論更合適,理論上,如果利用非線性波的勢函數(shù),建立非線性波作用下的滲流場解析解,文中方法可以推廣到非線性波情形。2)滲流計(jì)算中,沒有考慮海床的變型,這就意謂著模型只能用于沙質(zhì)海床。

基于Navier-Stokes方程,采用差分法或有限元法計(jì)算的波浪沖刷暴發(fā)時(shí),都用振蕩流代替波浪運(yùn)動(dòng),可能與真實(shí)的波浪運(yùn)動(dòng)有差異。而邊界元方法則可以準(zhǔn)確模擬波浪過程。對于波浪作用下的沖刷暴發(fā),邊界元方法可能更為簡單而準(zhǔn)確。當(dāng)然,邊界元方法的局限性也是明顯的,對于單向流或波流聯(lián)合作用下的沖刷暴發(fā),水體運(yùn)動(dòng)都不能用勢流模型模擬,邊界元也就失去了意義。

[1]LI F,CHENG L.Numerical model for local scour under offshore pipelines[J].Journal of Hydraulic Engineering,1999,125(4):400-406.

[2]LI F,CHENG L.Prediction of lee-wake scouring of pipelines in currents[J].Journal of Waterway,Port,Coastal,and Ocean Engineering,2001,127(2):106-112.

[3]BRORS B.Numerical modeling of flow and scour at pipelines[J].Journal of Hydraulic Engineering,1999,125(5):511-523.

[4]SUMER B M,JENSEN H R,FREDSOE 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.

[5]YEGANEH-BAKHTIARY A,KAZEMINEZHAD M H,ETEMAD-SHAHIDI A,et al.Euler-Euler two-phase flow simulation of tunnel erosion beneath marine pipelines[J].Applied Ocean Research,2011,33(2):137-146.

[6]DUPUIS A,CHOPARD B.Lattice gas modeling of scour formation under submarine pipelines[J].Journal of Computational Physics, 2002,178(1):161-174.

[7]MAO Y.The interaction between a pipeline and an erodible bed[D].Lyngby,Denmark:Technical University of Denmark,1986.

[8]CHIEW Y M.Mechanics of local scour around submarine pipelines[J].Journal of Hydraulic Engineering,1990,116(4):515-529.

[9]SUMER B M,FREDSOE J.Onset of scour below a pipeline exposed to waves[J].International Journal of Offshore and Polar Engineering,1991,1(3):189-194.

[10]KLOMP W H G,HANSEN E A,CHEN Z,et al.Pipeline seabed interaction,free span development[C].Proceeding 5th International Offshore and Polar Engineering Conference Netherlands:The Hague,1995(II):117-122.

[11]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.

[12]YANG B,GAO F P,WU Y X,et al.Numerical study of the occurrence of pipeline spanning under the influence of steady current[J]. Shipbuilding of China,2005,46(Suppl.1):221-226.楊兵,高福平,吳應(yīng)湘,等.海流引起海底管道懸空的數(shù)值模擬[J].中國造船,2005,46 (增刊1):221-226.

[13]LIANG D,CHENG L.A numerical model of onset of scour below offshore pipelines subject to steady currents[J].Frontiers in Offshore Geotechnics,2005:637-644.

[14]ZANG Z,CHENG L,ZHAO M,et al.A numerical model for onset of scour below offshore pipelines[J].Coastal Engineering,2009, 56:458-466.

[15]ANG W T.A Beginner’s course in boundary element methods[M].Boca Raton USA:Universal Publishers,2007:1-34

Boundary Element Method for Calculating the Onset of Scour Beneath Marine Pipelines Under the Wave Action

XIA Hua-yong1,WANG Jun-qin2
(1.South China Sea Marine Forecasting Center,Guangzhou 510300,China; (2.Engineering Design Department,CNOOC Research Institue,Beijing 100027,China)

The critical condition of the onset of scour beneath marine pipelines is one of the important fundamental parameters for the design and operation of submarine pipelines.Sumer et al gave the test curves of pipeline scour onset under the wave action,which included only three curves of Keulegan-Carpenter (KC)number.For an arbitrary KC number,however,the Sumer et al's method is not convenient for use and also not accurate.It is,therefore,necessary to establish a numerical model for calculating the critical condition of the onset of pipeline scour under the wave action.Because the intersection of the pipelines with the seabed can form a complex geometrical shape and the cross points are numerical singularity,the calculation of the critical condition of the onset of pipeline scour becomes very difficult.For this reason,the Boundary Element Method(BEM)is applied to calculate the wave motion and the seepage pressure.The boundary elements can precisely fit the boundary between the seabed and the pipeline wall,easily solve the problem about the numerical singularity at the cross points of the pipeline and the seabed and directly calculate the vertical pressure gradient at the cross points.Thus,the critical condition of the onset of pipeline scour can be accurately calculated.The calculated results agree well with the experiments available in literatures.

marine pipeline;onset of scour;boundary element method;numerical singularity

September 13,2016

P75

A

1002-3682(2017)02-0009-08

10.3969/j.issn.1002-3682.2017.02.002

2016-09-13

水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目——沙脊沙波及路由海床穩(wěn)定性研究(HESS1401)

夏華永(1967-),男,湖南沅江人,研究員,博士,主要從事海洋動(dòng)力學(xué)調(diào)查方面研究.E-mail:xiahuayong2001@21cn.com