夏華永,王俊勤

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

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

夏華永1,王俊勤2

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

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

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

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

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

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

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

(李 燕 編輯)

1 沖刷暴發(fā)機制

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

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

Zang等[14]取滲流出流點A處的垂向壓力梯度與管道下方的平均壓強相關

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

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

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

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

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

2 控制方程

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

對于微幅波,波浪勢為

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

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

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

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

3 邊界元法

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

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

圖2 波浪、滲流計算區(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)給定壓強,如果存在不滲水層,則取?P/?n=0。對于Laplace方程,其基礎解為

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

得到Laplace方程的解

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

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

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

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

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

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

2)用邊界元法計算t=iΔt,(i+1)Δt時刻的波浪勢函數(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,計算t=(i+1/2)Δt時刻海底的波浪壓強,并作為滲流的邊界條件,在邊界L2,L4上,根據(jù)式(9)或(10)給定?P/?n,在邊界L1上,如透水,則據(jù)式(9)給定壓強,如果不滲水,則取?P/?n=0;

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

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

圖5 Test O16中交叉點壓強梯度過程線[11]Fig.5 Time series of pressure gradient at the cross point in Test O16[11]

4 模型驗證

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

圖6 模擬交叉點A處壓強梯度過程Fig.6 Simulation of the time series of pressure gradient at cross point A

圖7 波浪質點速度與壓強水頭分布Fig.7 Distributions of the wave particle velocity and the pressure head

5 結 語

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

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

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

[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.楊兵,高福平,吳應湘,等.海流引起海底管道懸空的數(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

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

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