成小飛，王永學，王國玉，任 冰

（大連理工大學海岸和近海工程國家重點實驗室，大連 116024）

海底子母管線在規(guī)則波作用下的數(shù)值研究

成小飛，王永學，王國玉，任 冰

（大連理工大學海岸和近海工程國家重點實驗室，大連 116024）

采用三步有限元法離散N-S方程，建立了數(shù)值波浪水槽模型。數(shù)值研究了規(guī)則波作用下海底子母管結構的水動力特性。對海底子母管線所受波浪力的數(shù)值結果與物模實驗結果進行比較。基于該數(shù)值模型，還考察了不同的子母管間相對縫隙G/D（0.1，0.25，0.5，0.75）對海底子母管線水動力特性的影響。分析得到了不同G/D下海底子母管線的渦脫落特性以及管線水動力系數(shù)隨G/D的變化規(guī)律。結果表明，數(shù)值結果和實驗結果吻合較好，該數(shù)值波浪水槽模型可用來計算波浪作用在海底子母管線上的水動力。關于不同G/D，得到了2種不同的渦脫落模式，在G/D較小（0.1，0.25）時，為“反相同步脫落”模式，在G/D較大（0.5，0.75）時，為“同相同步脫落”模式。

海底子母管線；規(guī)則波；三步有限元；CLEAR-VOF；渦脫落

Biography：CHENG Xiao-fei（1985-），male，doctor student.

在實際工程中，不等直徑的海洋輸油氣管線有時因為技術和經濟的原因，經常以一定間隙捆在一起形成管束，置于海底。管束一般由一根母（大）管和幾根子（?。┕芙M成。一種比較流行的管束形式是海底子母管結構［1］，它由一根母管和一根子管組成，子管以一定間隙置于母管正上方。海底子母管結構與單管相比，其周圍的流場形態(tài)及水動力特征更加復雜，所以關注這種海底子母管結構的水動力特性至關重要。

關于海底子母管線這種特殊結構，國內外很多學者已做了一些研究，主要集中在穩(wěn)定流中，比如Kalghatgi等［2］采用拖車實驗的方法研究得到了穩(wěn)定流情況下，管線與海床的相對間隙、Re（Reynolds）數(shù)的變化對海底子母管線水動力系數(shù)的影響；Kamarudin等［3］、Zhao等［4］數(shù)值研究了穩(wěn)定流情況下海底子母管線水動力問題，研究得到了不同Re數(shù)下子母管周圍的壓力分布，水動力系數(shù)以及渦脫落形式。

目前關于海底子母管線在波浪作用下的研究成果主要來自物模實驗，比如馬良等［5］和李玉成等［6］在波流水槽內研究了波流共同作用下KC（Keulegan-Carpenter）數(shù)對海底子母管線水動力系數(shù)的影響；Brankovic等［7］通過拖車造波的方式研究了波流共存下海底子母管線的水動力；Cheng等［8］在波浪水槽內著重研究了規(guī)則波和不規(guī)則波作用下海床對子母管線水動力的近壁影響。而目前仍缺乏關于海底子母管線在波浪作用下的數(shù)值研究。

本文采用三步有限元方法來離散二維粘性不可壓流體的雷諾平均N-S方程和連續(xù)性方程，用k-ω模型來模擬流體的湍流運動，采用CLEAR-VOF方法來追蹤運動水體的自由表面，構建了二維數(shù)值波浪水槽模型［9-10］。用該數(shù)值波浪水槽進行了海底子母管線在規(guī)則波浪場中受力的數(shù)值實驗，并同物模試驗結果進行比較驗證，著重考察了子、母管間相對縫隙G/D（0.1，0.25，0.5，0.75）對海底子母管線水動力特性的影響。

1 基本控制方程

基本控制方程是二維粘性不可壓雷諾平均N-S方程和連續(xù)性方程。在笛卡爾坐標系中無量綱化后的基本方程如下

式中：xi（i=1，2）為水平或垂直方向坐標；ui為i方向的速度分量；p為壓力；k為湍動能；Re（=U0D/ν）為雷諾數(shù)，由母管中心位置處水質點的最大水平向速度U0、母管直徑D和流體運動粘性系數(shù)v計算得到；t為無量綱化時間；νt為渦粘性；Sij為平均流動變形率張量，其方程如下

k-ω湍流模型［11-12］用來計算湍動能k和渦粘性νt，其方程如下

式中：ω=ε/（β*k）為湍動能比耗散率，ε為湍動能耗散率；pk=2νtSij?ui/?xj為湍動能生成項；渦粘性νt由下式計算得到

式（4）～式（6）中各模型參數(shù)定義如下

2 數(shù)值計算方法

本文對基本控制方程的求解采用三步有限元法［13］，即將時間導數(shù)按泰勒級數(shù)展開，一個時間步內分3次進行迭代，以獲得穩(wěn)定的數(shù)值解。與SUPG等傳統(tǒng)迎風有限元方法相比，三步有限元法數(shù)值解精度高（三階），數(shù)值實現(xiàn)過程簡單，并具有更高的計算效率。N-S方程（1）應用三步有限元法，可得到如下時間離散形式

式中：上標 n，n+1/3，n+1/2 和 n+1 分別表示 n，n+Δt/3，n+Δt/2 和 n+Δt時刻的量值。式（7）～式（9）可用標準的Galerkin加權余量有限元法進行空間導數(shù)離散。

針對本文的湍流動問題，式（7）～式（9）中k和vt不是常數(shù)，需要進行線性化處理，即利用上一時刻的流場計算結果，根據(jù)k-ω湍流模型式（4）～式（6）求得k和vt，并在當前時間步內將之視為常數(shù)添加到式（7）～式（9）中。其中對于k-ω湍流模型，其控制方程也是對流擴散方程類型，形態(tài)上同于N-S方程，因此對式（4）～式（6）的離散和求解也可參照N-S方程的處理方法。

本文建立數(shù)值波浪水槽，可反映出計算域中自由表面的實時變化特征。在任意時間步內，通過三步有限元法數(shù)值求解N-S方程，得到該時間步對應的速度場和壓力場后，對計算域內的自由表面進行重構。關于自由表面的重構本文采用的是Ashgriz等［14］提出的CLEAR-VOF（Computational Lagrangian-Eulerian Advection Remap Volume of Fluid）方法。它是以有限元為基礎的一種高效的自由表面重構技術。CLEAR-VOF方法是在歐拉網格基礎上采用了拉格朗日思想追蹤單元內水團在新時刻的位置以求解每個單元內VOF函數(shù)。那么在經歷一個時間步Δt后，水團頂點的位移和新位置可表示為

水團幾何頂點新位置確定后，水團位置也由此得到。這種方法的主要優(yōu)點包括無需求解關于流體體積函數(shù)的微分方程，而是在拉格朗日意義下實現(xiàn)流體的對流輸運；界面重構對具體的網格劃分方式無特殊要求，不但適用于規(guī)則網格，同時也可用于非規(guī)則網格。

3 數(shù)值波浪水槽

圖1 數(shù)值波浪水槽計算域Fig.1 Computational domain of numerical wave flume

數(shù)值波浪水槽計算域如圖1所示，水槽長6 m，高0.4 m，水深0.3 m。海底子母管結構如圖2所示，位于水槽中央，母管直徑D為0.02 m，子管直徑d為0.008 m，母管與海底之間的間隙比e/D=0.25，子母管之間的相對縫隙G/D=0.25，海底及管線表面均光滑，不考慮粗糙率的影響。

對計算域采用四邊形線性等參元進行網格劃分。圖3給出了子母管結構周圍的網格劃分詳細圖。整個計算域網格共有61 530個節(jié)點，60 784個單元，網格的最小尺寸約為0.000 5 m，最大尺寸約為0.03 m，時間步長為0.001 s。

計算域左邊界和右邊界分別設為可吸收式無反射造波邊界與開邊界［15］，海底和管線表面為無滑移邊界條件，水平和垂向速度均為零，計算域頂表面為滑移邊界條件，垂向速度為零。初始時刻，流體是靜止的，整個計算域的速度和壓力都設為零。

圖2 海底子母管結構Fig.2 Submarine piggyback pipeline configuration

圖3 子母管周圍的網格（e/D=G/D=0.25）Fig.3 Meshes near the piggyback pipeline（e/D=G/D=0.25）

4 數(shù)值計算結果和實驗結果的比較

為了驗證數(shù)值波浪水槽模型的可靠性，本文將數(shù)值計算的結果和物模實驗結果進行比較。物模實驗在大連理工大學海岸與近海工程國家重點實驗室海洋環(huán)境水槽中進行，通過水下雙向測力傳感器同步測得了海底子母管系統(tǒng)在規(guī)則波作用下總的水平力和升力（垂向力）。關于物模實驗的細節(jié)另文給出。為了便于比較，物模實驗工況和數(shù)值計算的工況選為一致。

圖 4給出了波高 H=0.1 m，周期 T=1.7 s，水深 h=0.3 m，e/D=G/D=0.25條件下（KC=25.14，Re=5 880），由數(shù)值波浪水槽計算得到的波面歷時過程線（波面過程線的采樣點在管線中心位置正上方）、子母管系統(tǒng)所受總的水平力和升力歷時過程線與實測結果的比較。由圖4可知，數(shù)值計算得到的波面過程線、水平力過程線與實測數(shù)據(jù)結果吻合較好，兩者在幅值上基本相同，相位也基本吻合。兩者升力過程線在幅值和相位上也大體一致，但是在一個波周期內，實測升力過程線比數(shù)值計算升力過程線多一個或幾個小的峰（谷）值，這是由于實驗條件下升力的隨機性很大，波浪場中微小的擾動即可造成升力的變化。從整段升力過程線來看，數(shù)值計算與實測結果基本吻合。綜上分析，本文采用的數(shù)值計算方法是可靠的，建立的數(shù)值波浪水槽模型可以用來計算波浪作用下海底子母管線的受力。

圖4 數(shù)值計算結果與物模實驗結果的比較（KC=25.14，e/D=G/D=0.25）Fig.4 Comparison of numerical and experimental results（KC=25.14，e/D=G/D=0.25）

5 G/D對海底子母管線水動力特性的影響

基于本文建立的數(shù)值波浪水槽模型，考察了子母管間的相對縫隙G/D（0.1，0.25，0.5，0.75）對海底子母管線水動力特性的影響，分析得到了在不同G/D下海底子母管線的渦脫落特性以及管線水動力系數(shù)隨G/D的變化規(guī)律。在數(shù)值計算中，取子管與母管管徑之比d/D為0.4，母管與海床間隙比e/D為0.25，波浪為規(guī)則波（波高 H=0.1 m，周期 T=1.7 s），相應 KC=25.14。

5.1 渦脫落特性

圖5分別給出了不同G/D（0.1，0.25，0.5，0.75）下，海底子母管線（子管、母管、子母管系統(tǒng)）在一個波浪振蕩周期內（t/T=7.4～8.4）所受水平力Fx和升力FL的歷時過程線。從圖5可以看出，子管、母管和子母管系統(tǒng)所受水平力Fx在4種G/D下均同步振蕩，與波浪振蕩頻率相同，但是所受升力FL比較復雜，其振蕩頻率是兩倍甚至三倍的波浪振蕩頻率。在G/D較小時（G/D=0.1），子管和母管所受升力相互排斥，反相振蕩，這是因為G/D較小時子、母管之間的相互干擾較明顯。隨著G/D的增大，子、母管間的相互干擾逐漸減弱，子管和母管升力的振蕩趨于同步。子母管系統(tǒng)所受波浪力的歷時過程線與母管的非常接近，這主要跟母管管徑相對較大，在子母管系統(tǒng)中占主導地位有關。

圖5 不同G/D下海底子母管線水平力和升力的歷時過程線Fig.5 Time histories of in-line and lift forces on submarine piggyback pipeline for different G/D

在圖5中不同G/D（0.1，0.25，0.5，0.75）對應的波浪力歷時曲線，取其前半周期升力曲線出現(xiàn)的2個極值時刻 a、b，考察其管線周圍渦量的分布特性。圖 6 分別給出了不同 G/D（0.1，0.25，0.5，0.75）下，在瞬時 a、b（對應于圖 5 中的瞬時 a、b）海底子母管線周圍渦量的分布圖。渦量 wz由公式 wz=0.5×（?u2/?x1-?u1/?x2）計算得到，圖6中的虛線表示渦量值為負（負渦），實線表示渦量值為正（正渦）。由圖6-a、6-b的渦量分布圖可見，子母管間縫隙較小，子管與母管縫隙間的渦因受到抑制較小，而子管頂部和母管底部的渦不斷增大，在時刻a，子管頂部的負渦和母管底部的正渦同時在管線右側脫落，引起子、母管升力達到極值，而且兩者受力是反相的（如圖5-a、圖5-b），稱為“反相同步脫落”模式。隨著波浪水質點速度的轉向，母管右側的正渦先轉向反作用在母管上，引起母管升力達到極值。隨后在時刻b，子管右側脫落的負渦也發(fā)生轉向，同時作用在子管和母管上，引起子、母管所受水平力和升力均達到極值。

圖6 不同G/D下海底子母管線周圍渦量的分布圖Fig.6 Vorticity contours around submarine piggyback pipeline for different G/D

由圖6-c、圖6-d的渦量分布圖可見，在G/D較大時（0.5，0.75），子管與母管縫隙間的渦所受抑制逐漸減弱，子母管間縫隙處的渦逐漸增大。在時刻a，子管頂部的負渦和母管頂部的負渦均在管線右側同步脫落，引起子、母管所受升力達到極值，并且由圖5-c、圖5-d可以看出，子、母管所受升力同相，稱為“同相同步脫落”模式。隨后在時刻b，由于波浪水質點速度的轉向，子管右側脫落的負渦和母管右側脫落的負渦均發(fā)生轉向，分別反作用在各自的管線上，引起各自管線所受升力達到極值。

由以上分析，可以得到波浪作用下，管線周圍渦的脫落或者脫落渦的反作用都會引起管線所受水平力或者升力達到極值。在G/D=0.1，0.25時，子管對母管的影響較大，母管不僅受到母管自身脫落渦的反作用產生升力極值，同時還受到子管脫落渦的反作用產生極值。而隨著G/D的增大，子管對母管的影響逐漸減小，在G/D=0.5，0.75時，母管僅受到自身脫落渦的反作用。因此從渦的脫落特性，可以說明圖5中在G/D=0.1，0.25時，母管所受升力的峰（谷）值個數(shù)要多于G/D=0.5，0.75時的情形。

圖 7 分別顯示了不同 G/D（0.1，0.25，0.5，0.75）對應的子管、母管和子母管系統(tǒng)的無量綱Strouhal數(shù)。St=f0D0/U0，其中f0為管線升力的振蕩頻率；D0為管線直徑，子、母管和子母系統(tǒng)分別取d，D和等效直徑DE（DE=D+d+G）；U0為來流速度，本文取母管中心位置處波浪的水質點水平速度。從圖7可以看出，在G/D=0.1，0.25時，子管僅存在單一St數(shù)，而母管和子母管系統(tǒng)存在2個St數(shù)，這是因為在G/D較小時，子管對母管的影響較大，子管脫落的渦對母管也有額外的反作用。隨著G/D的增大，子管脫落的渦對母管的反作用減弱，在G/D=0.5，0.75時，子管、母管和子母管系統(tǒng)均只存在單一St數(shù)。另外，由圖7中各管線St數(shù)的比值關系還可以發(fā)現(xiàn)，子、母管和子母管系統(tǒng)升力振蕩的主頻率是相同的，均是波浪振蕩頻率的2倍，而在G/D較小時還存在3倍于波浪振蕩頻率的次頻。

圖7 不同G/D下海底子母管線的St數(shù)Fig.7 Stnumbers of submarine piggyback pipeline for different G/D

5.2 水動力系數(shù)

在波浪作用下，管線所受水平力和升力可按Morison方程計算得到

式中：Fx（t）為管線所受的水平波浪力；FLmax為管線所受的最大升力；ρ為流體密度；l為管線的長度；u（t）、a（t）分別為流場未受擾動時母管中心位置處波浪的水質點水平速度和加速度；umax為流場未受擾動時母管中心位置處波浪的最大水質點水平速度；D0為管線直徑，子、母管和子母系統(tǒng)分別取d，D和等效直徑D（EDE=D+d+G）；CD、CM和CL分別為管線的拖曳力系數(shù)、慣性力系數(shù)和升力系數(shù)。

CD和CM可由最小二乘法計算得到，而升力系數(shù)C（L，升力方向背離海床；，升力方向指向海床）由極值法計算得到。關于水動力系數(shù)的詳細分析方法可參考李玉成等［16］。圖8分別給出了在規(guī)則波作用下（KC=25.14）海底子母管結構中子管、母管和子母管系統(tǒng)的水動力系數(shù)CD，CM，CL隨間隙比G/D的變化規(guī)律。

從圖8可以看出，子管、母管和子母管系統(tǒng)的CD，CM關于G/D的變化規(guī)律大體一致，均隨G/D的增大而減小，并且母管和子母管系統(tǒng)的CD，CM值比較接近，這跟母管在子母管系統(tǒng)中占主導地位有關。子管的CD，CM值要小于母管和子母管系統(tǒng)的CD，CM值，是后者的65%～90%。子管和母管的升力系數(shù)CL關于G/D的變化規(guī)律是相反的，子管的隨G/D的增大而減小，的絕對值隨G/D的增大而增大，這是因為在G/D較小時，母管的存在抑制了子管底部渦的發(fā)展，導致子管很小較大，隨著G/D的增大，母管的影響逐漸減弱絕對值不斷增大，不斷減小。相反，母管的隨G/D的增大而增大的絕對值隨G/D的增大而減小。另外，子母管系統(tǒng)的升力系數(shù)隨G/D的變化趨勢同母管大體一致，量值上也比較接近。

圖8 海底子母管線水動力系數(shù)隨G/D的變化規(guī)律Fig.8 Variation of hydrodynamic coefficients on submarine piggyback pipeline with G/D

6 結論

本文采用三步有限元法來離散N-S方程，k-ω模型來模擬流體的湍流運動，CLEAR-VOF方法來追蹤運動水體的自由表面，建立了規(guī)則波作用在海底子母管線上的數(shù)值波浪水槽模型，將數(shù)值計算結果同物模實驗結果進行了比較驗證。并考察了在KC=25.14條件下，子母管間的相對縫隙G/D（0.1，0.25，0.5，0.75）對海底子母管線水動力特性的影響。其主要結論如下：

（1）經比較驗證，本文所建立的數(shù)值波浪水槽模型可以用來計算海底子母管線在規(guī)則波作用下的水動力。

（2）關于不同G/D，本文得到了海底子母管線在規(guī)則波作用下2種不同的渦脫落模式。在G/D=0.1，0.25時，存在一種渦脫落模式稱為“反相同步脫落”模式；在G/D=0.5，0.75時，存在另一種渦脫落模式稱為“同相同步脫落”模式。

（3）子管、母管和子母管系統(tǒng)的C，C關于G/D的變化趨勢大體一致，均隨G/D的增大而減小。子管的DM隨G/D的增大而減小，的絕對值隨G/D的增大而增大，母管的升力系數(shù)隨G/D的變化趨勢與子母管系統(tǒng)相一致，但與子管的變化趨勢相反。

［1］楊琥，倪浩，朱曉環(huán).一種新型的置換海底子母管道技術［J］.中國造船，2007（48）：563-570.

YANG H，NI H，ZHU X H.An applicable replacement bundled pipeline structure for offshore marginal oilfield development［J］.Shipbuilding of China，2007（48）：563-570.

［2］Kalghatgi S G，Sayer P G.Hydrodynamic forces on piggyback pipeline configurations［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，1997，123（1）：16-22.

［3］Kamarudin M H，Thiagarajan K P，Czajko A.Analysis of current-induced forces on offshore pipeline bundles［C］//Witt P J.Proceeding of 5th International Conference on CFD in the Process Industries.Melbourne：CSRIO Australia，2006.

［4］ZHAO M，CHENG L，TENG B.Numerical modeling of flow and hydrodynamic forces around a piggyback pipeline near the seabed［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，2007，133（4）：286-295.

［5］馬良，王金英，孫紹述，等.海底（子-母）管道在波浪作用下水動力載荷的實驗研究［J］.大連大學學報，1993，3（4）：54-63.

MA L，WANG J Y，SUN S S，et al.Experimental study of hydrodynamic force loading under the action of wave current on submarine composite pipelines［J］.Journal of Dalian University，1993，3（4）：54-63.

［6］李玉成，張寧川，孫姎.波流共同作用下近底子母管線的水動力特征［J］.水動力學研究與進展，1994，9（1）：51-59.

LI Y C，ZHANG N C，SUN Y.The hydrodynamic characteristic of submarine composite pipeline in wave-current coexisting field［J］.Journal of Hydrodynamics，1994，9（1）：51-59.

［8］CHENG X F，WANG Y X，WANG G Y.The effect of the seabed proximity on the hydrodynamic forces of the piggyback pipeline under wave action［C］//Chung J S.Proceeding of 30th International Conference on Ocean，Offshore and Arctic Engineering.New York：ASME，2011.

［9］LU L，LI Y C，TENG B，et al.Numerical simulation of turbulent free surface flow over obstruction［J］.Journal of Hydrodynamics，2008，20（4）：414-423.

［10］孫英偉，陳兵，康海貴.遠破波作用數(shù)值模擬的 CLEAR-VOF 模型［J］.水科學進展，2010，21（6）：795-800.

SUN Y W，CHEN B，KANG H G.Numerical simulation of broken wave with the CLEAR-VOF-FEM model［J］.Advances in Water Science，2010，21（6）：795-800.

［11］Wilcox D C.Reassessment of the scale-determining equation for advanced turbulence models［J］.AIAA Journal，1988，26（11）：

1 299-1 310.

［12］Wilcox D C.Simulation of transition with a two-equation turbulence model［J］.AIAA Journal，1994，32（2）：247-255.

［13］JIANG C B，Kawahara M.The analysis of unsteady incompressible flows by a three-step finite element method［J］.International Journal of Numerical Methods in Fluids，1993，6：793-811.

［14］Ashgriz N，Barbat T，WANG G.A computational Lagrangian-Eulerian advection remap for free surface flows［J］.International Journal of Numerical Methods in Fluids，2004，44：1-32.

［15］王永學.無反射造波數(shù)值波浪水槽［J］.水動力學研究與進展，1994，9（2）：205-213.

WANG Y X.Numerical wave channel with absorbing wave-maker［J］.Journal of Hydrodynamics，1994，9（2）：205-213.

［16］李玉成，陳兵，王革.波浪對海底管線作用的物理模型實驗及數(shù)值模擬研究［J］.海洋通報，1996，15（4）：58-65.

LI Y C，CHEN B，WANG G.Physical model test and numerical simulation of pipeline under wave action［J］.Marine Science Bulletin，1996，15（4）：58-65.

Numerical study of submarine piggyback pipeline under regular wave action

CHENG Xiao-fei，WANG Yong-xue，WANG Guo-yu，REN Bing
（State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian116024，China）

A numerical wave flume model was established，in which the Navier-Stokes equation was dispersed using the three-step finite element method.Then the hydrodynamic characteristics of the submarine piggyback configuration under regular wave action were numerically studied.The numerical results for wave forces on submarine piggyback pipeline were compared with the experimental results.Based on this numerical model，the effect of the spacing ratios between the small and the large pipelineG/D（0.1，0.25，0.5，0.75）on the hydrodynamic characteristics of the submarine piggyback pipeline was also investigated.The vortex shedding characteristics of the piggyback pipeline about differentG/Dand the variation of hydrodynamic coefficients of the pipeline withG/Dwere presented in this paper.The results show that both numerical and experimental results are in good agreement,and the numerical wave flume model can be applied to predict the hydrodynamic forces on the submarine piggyback pipeline under wave action.Two different vortex shedding patterns are found about differentG/D.One is the“anti-phase-synchronized”pattern forG/D=0.1，0.25，the other is the“in-phase-synchronized”pattern forG/D=0.5,0.75.

submarine piggyback pipeline；regular wave；three-step finite element;CLEAR-VOF；vortex shedding

TV 139.2；O 242.1

A

1005-8443（2012）03-0185-09

2012-01-05；

2012-03-23

國家自然科學創(chuàng)新研究群體基金（50921001）；國家重點基礎研究發(fā)展計劃（973計劃）（2011CB013702）

成小飛（1985-），男，江蘇省南通人，博士研究生，主要從事波浪與海洋結構物的作用研究。