基于勢流理論的內孤立波與立管作用數值研究1)

胡英杰 鄒 麗,?,2) 孫 哲 金國慶 馬鑫宇

* (大連理工大學船舶工程學院,遼寧大連 116024)

? (工業裝備結構分析國家重點實驗室,遼寧大連 116024)

引言

海洋內波是一種發生在海面以下的波浪,具有波幅大、周期長的特征,內波在海洋中經常以大波幅的內孤立波及波列的形式傳播,在世界范圍的海域內廣泛存在[1-3].在傳播過程中,內孤立波的流場會在海洋表面產生明顯的特征,因而使用合成孔徑雷達(SAR)對海面特征的捕捉成為對其觀測的重要依據[4-6].同時,使用聲學多普勒流速儀(ADCP)等儀器的海洋實測也在多個區域觀測到了內孤立波的存在[7-8].通過遙感和實地觀測表明我國南海也是內孤立波頻發的海域[9].在海洋資源開發過程中,海洋石油和深海礦產輸送立管長期在海上作業,海洋立管在內孤立波流場的作用下會產生劇烈而大幅度的動力響應,對管道的構形、結構強度、作業效率與安全都產生嚴重影響[10].海洋內孤立波不僅對水下潛器和結構物影響巨大,還會對海洋平臺等水面結構物產生影響.內孤立波在遠距離傳播過程中,會在海洋表面引導出明顯的海水流動,因而會造成海洋平臺的長距離漂移,對錨泊結構是一個重大考驗,影響正常生產作業,造成嚴重的經濟損失[11].在安德曼海作業的一個石油鉆井機曾被海洋內孤立波的流場推移了30.48 m,嚴重影響作業安全[12].

近年來,內孤立波對海洋立管的作用引起了更多學者的關注.文獻[13-15]針對內孤立波對海洋立管的作用問題開展了系統深入的研究,通過KdV 和mKdV 方程計算了內孤立波的流場特征,采用動力有限元模型得到了頂張力立管的動力響應,并討論了內波參數對管道響應的影響.文獻[16-18]采用試驗和數值模擬方法研究了內孤立波與海流共同作用下管道與海洋平臺樁柱的動力響應問題,獲得了不同入射角度及不同波浪參數條件下的管道變形特征.Wang等[19]采用非線性有限元模型計算了弱約束深海礦產輸送管道在內孤立波作用下的變形與運動特性,討論了內孤立波不同入射角度及水面平臺運動對管道的影響.Wang等[20]采用試驗方法研究了內孤立波對小尺度樁柱的作用力分布特征,討論了不同水深條件下樁柱的受力情況,分析了水動力系數選取與KC數的關系.Fan等[21]基于歐拉-伯努利梁理論,考慮了海流、表面波和內孤立波的共同作用,發現內孤立波的存在明顯增大了立管的變形包絡,并且分層流體的密度差也會對立管變形產生影響.李朝瑋等[22]采用準靜態分析法研究了非均勻海流、表面波和內孤立波作用下隔水管的變形,發現內孤立波對立管的橫向位移和應力都有不可忽視的影響且表面波和內波相位角分別在12°和0°時對立管的影響最大.Zhang等[23]采用數值模擬和試驗方法研究了平臺運動和內孤立波共同作用下頂張力立管的渦激振動特性,發現內孤立波的波幅越大、平臺升沉運動的頻率越高,立管的渦激振動越明顯.Wang等[24]在分層流水槽中開展模型試驗分析了斜向入射的內孤立波的波流結構特征及其對水下細長結構的載荷作用特性.文獻[25-28]采用莫里森經驗公式、模態分解和回歸分析等方法對內孤立波對柱體及管道的動力載荷進行了深入分析.

內孤立波流場速度和加速度信息的準確獲取是立管載荷確定的關鍵,進而影響整個計算模型的準確性.對于大尺度的海洋內孤立波,無法開展實驗獲得流場信息,同時求解NS 方程的CFD 方法計算量較大.目前對大尺度內孤立波與立管作用的研究中,內孤立波流場信息多是基于解析模型近似給定的.本文工作基于勢流理論采用多域邊界單元法建立數值模型求解了內孤立波的流場特征,較為高效又能實時準確地計算流場信息.依據流場信息采用莫里森方程給定立管載荷,采用動力學有限元方法計算了頂張力立管和弱約束立管的動力響應,以期為海洋立管設計在考慮內孤立波作用時提供參考.

1 計算模型

1.1 分層流體運動勢流計算模型

內孤立波與立管作用如圖1 所示,兩層流體的深度分別為h1,h2;密度分別為ρ1,ρ2.下凹型內孤立波位于密度界面分層處.

圖1 內孤立波與立管作用示意圖Fig.1 Schematic diagram of the interaction between internal solitary wave and the riser

兩層流體都假定為無旋,無黏性,不可壓縮的,因此兩層流體中存在流動勢函數φk滿足拉普拉斯方程

在分層流體界面處的運動學邊界條件為上下層流體沿界面法線方向的速度連續,因此有

其中η2為內孤立波波面,n為計算區域Ω1,Ω2在界面處的外法線方向矢量.

在分層流體界面處的動力學邊界條件為上下層流體區域在界面處的壓力連續,因此有

其中g 為重力加速度.海底滿足不可穿透條件

由于內孤立波傳播時在海洋表面產生的波動較小,為了提高計算效率,本文對上表面采用剛蓋假定,因此上表面的邊界條件為

通過格林第二定理,可將上下兩層流體區域Ω1,Ω2的控制方程拉普拉斯方程轉化為兩個邊界上的曲線積分方程

其中r為源點,q為場點,G(r,q)=-1/(2π)ln|r-q|為格林函數,c(r)=1-α(r)/(2π),α(r)為場點r處的單元夾角,Γk為上下層區域邊界.

由于內孤立波位于密度分層的流場中,所以計算過程需要涉及兩個流體區域,進而要將單獨區域的邊界單元法擴展到多個區域[29-30].對于計算區域Ω1,Ω2的邊界分別采用n個線單元進行離散,從而積分方程 (6) 可以寫成如下的形式

根據兩個區域的關系,可以將邊界點分為公共邊界點和自身邊界點,從而將積分方程進一步處理為

其中,上標c和s分別指公共邊界點和自身邊界點.

波面位置和速度勢的更新格式分別根據運動學和動力學邊界條件給定

內孤立波初始波形的給定依據強非線性的MCC 模型,實現對大波幅內孤立波的準確模擬.

1.2 管道動力響應模型

在內孤立波由遠方向立管傳播過程中,立管處的流場速度及加速度是實時變化的,因此不能用簡單的均勻流或剪切流靜力假定求解,需要采用動力學進行求解.

海洋立管的動力學控制方程

其中,M,C,K分別為質量陣、剛度陣和阻尼陣,x(t),分別為管道的變形位移、速度和加速度.Q(t)為載荷向量,根據內孤立波流場的速度、加速度以由莫里森方程給定.

在計算中需要考慮管道軸向張力的影響,因此管道的剛度陣由兩個部分組成,分別是彈性剛度陣KE和幾何剛度陣KG,即

其中L為管道單元的長度.

阻尼陣采用瑞利阻尼,質量陣與剛度陣的線性組合形式

其中 α 和 β 為比例系數,根據模態頻率給定,分別為

ω1和ω2分別為管道的一階和二階頻率,ζ為管道的阻尼比.

采用Newmark-β法求解動力學控制方程 (11),獲得管道的動力響應.

2 內孤立波載荷

在通過對邊界積分方程的求解獲得區域邊界上的速度勢及其法向導數值后,流場內部速度勢可以用方程 (15) 表達.對積分方程 (15) 進行微分,即可將速度勢的積分方程轉化為關于流場速度的積分方程 (16),對該方程進行離散求解即可獲得內孤立波流場內部的速度分布特征,即

在進行計算求解時,當內部的待求點與邊界節點接近時,格林函數會表現出近奇異特征,導致計算精度降低甚至使計算失效.本文的研究中采用了單元子分法[31-32]來解決這一問題,得到了良好的計算結果.為了驗證計算模型的準確性,與文獻中的內孤立波水平速度垂向分布試驗結果[33]進行了對比.試驗在分層流體水槽中開展,上層流體深度h1=15 cm,下層深度h2=62 cm,密度分別為ρ1=0.999 kg/m3,ρ2=1022 kg/m3,內孤立波無量綱波幅為a/h1=0.22.計算結果與分層流體內孤立波傳播試驗流場的測量結果吻合良好如圖2 所示.計算得到的內孤立波流場速度分布如圖3 所示,計算結果表明內孤立波流場速度具有明顯的剪切特征,上層流體的質點速度與內孤立波的傳播方向相同,而下層流體則與之相反,且上層流體質點速度的絕對值要明顯大于下層流體.圖4 為內孤立波流場的水平速度在上層流體區域z/h1=0.2 處沿水平方向的分布情況.可以發現水平速度的最大值位于波谷位置(x/h1=25),且隨著與波谷距離的增大而迅速衰減.圖5為z/h1=0.2 處流體質點加速度在波長方向的分布情況,加速度的最大值位于波谷兩側,波谷位置處的加速度為零.

圖2 水平速度垂向分布與文獻[33]結果對比Fig.2 Comparison between the calculation result about the vertical distribution of horizontal velocity with the result of Ref.[33]

圖3 內孤立波流場速度分布Fig.3 Velocity distribution of internal solitary wave flow field

圖4 內孤立波水平速度沿水平方向的分布Fig.4 Horizontal velocity distribution of internal solitary wave

圖5 內孤立波水平加速度分布Fig.5 Horizontal acceleration distribution of internal solitary wave

由于深海礦產輸送管道具有大長細比特征,因此采用經典的莫里森公式作為橋梁將非線性內孤立波數值計算模型和管道響應的計算模型連接起來,從而實現內孤立波對管道作用的準確求解.

由于內孤立波對管道的作用過程中管道變形較大,所以要考慮管道的位移以及管道傾角的影響,因此內孤立波對管道作用力的莫里森公式寫成如下形式

其中,Cm為慣性力系數,本文計算取為2.0;Cd為拖曳力系數,取為1.2;Un為垂直于軸向的內孤立波流場速度,為對應的加速度;ρ為海水密度,x˙ 和x¨ 分別為管道運動的速度和加速度.根據前面內孤立波流場速度和加速度的計算結果及莫里森方程可知,當管道的直徑相對較小時,對管道響應起主要作用的載荷為內孤立波對管道的拖曳力,慣性力的影響相對較小.

3 管道動力響應

3.1 頂張力立管動力響應

在海洋石油開采等領域,頂張力立管得到了廣泛的應用,本文計算了內孤立波對頂張力立管的作用特性.為了驗證本文計算模型的準確性,將計算結果與文獻[34]結果進行了對比.上下層流體深度分分別為60 m 和412 m,密度分別為ρ1=1025 kg/m3,ρ2=1028 kg/m3,管道內外徑分別為0.26 m 和0.21 m,內孤立波波幅為75 m.對比結果如圖6 所示,可以發現兩者吻合較好.由于本文內孤立波流場的給定方法與文獻有一定區別,所以動力響應結果稍有偏差,管道順流向最大位移處的偏差約為3.34%.

圖6 管道流向最大位移與文獻[34]計算結果對比Fig.6 Comparison of the maximum displacement of the pipe in flow direction with the calculation results in the Ref.[34]

本文計算的頂張力立管的幾何參數和物理參數見表1.為了討論頂張力大小對管道動力響應的影響,頂張力T分別取為400 kN,500 kN,600 kN,700 kN,800 kN.

表1 頂張力立管參數Table 1 Parameters of the top tension riser

計算中采用的分層流體參數見表2.為了研究內孤立波參數對管道動力響應的影響,內孤立波的波幅分別給定為40 m,50 m,60 m,70 m,80 m.

表2 流體分層參數Table 2 Parameters of stratified fluids

不同內孤立波波幅時管道的變形特征如圖7 所示,可以發現立管的動力響應對波高敏感,隨著波高的增加,管道的流向位移明顯增大.同時,由于上層流體速度明顯大于下層,且在所研究問題中拖曳力遠大于慣性力,因此管道順流向的最大位移發生在上層區域.

圖7 不同內孤立波入射波幅時頂張力立管的流向變形特征Fig.7 Characteristics deformation of top tension riser in flow direction with different incident amplitudes of internal solitary waves

當內孤立波波幅為60 m 時不同頂張力條件下的立管變形特征如圖8 所示,可以發現頂張力對管道的流向變形具有重要影響.在相同內孤立波參數的條件下,施加的頂張力越大,管道的變形越小.圖9展示了不同張力時管道流向最大變形的變化趨勢.這是由于隨著頂張力T的增加,作用在管道單元上的有效張力Te也逐漸增大,通過方程 (13) 可知,管道的幾何剛度陣KG值將隨之增大,進而提高管道的整體剛度.因此,在相同參數的內孤立波作用下,頂張力越大,管道的順流向變形越小.

圖8 不同頂張力時管道的變形特征Fig.8 Deformation characteristics of the riser under different top tensions

圖9 管道最大流向位移隨頂張力大小的變化關系Fig.9 The relationship between the maximum displacement of the riser in flow direction and the value of the top tension

圖10 展示了z=-40 m 位置處管道的應力在不同內孤立波波幅條件下的變化特征.管道應力隨著波幅的增大而顯著增大,同時波幅較小的內孤立波作用時管道的最大應力時刻出現滯后.根據分層流體K d V 理論,內孤立波的波速可以表示為,其中c0和c1是與分層流體厚度和密度有關的量,因此當分層參數確定之后,內孤立波的波速與波幅密切相關,波幅越大的內孤立波具有更大的波速.因此,當波幅較小時,內孤立波傳播到管道位置處的時間顯著增加,管道的最大應力出現的時刻明顯滯后.

圖10 不同波幅內孤立波作用下管道的應力變化Fig.10 Stress changes of the riser under the action of solitary waves with different amplitudes

3.2 弱約束立管動力響應

本文討論的弱約束立管是指上端與水面設備剛固連接,下端為自由端的海洋立管,一般在深海采礦等領域得到應用.計算所采用的管道幾何和物理參數如表3 所示.

表3 弱約束立管參數Table 3 Parameters of the weakly constrained riser

為了討論內孤立波參數對弱約束立管的作用規律,計算了不同波幅內孤立波作用下的立管動力響應如圖11 所示.從圖中可以發現內孤立波波幅的變化對立管變形作用明顯,立管的流向位移隨初始內孤立波波幅的增大而顯著增大.

圖11 不同波幅內孤立波作用下弱約束立管的變形特征Fig.11 Characteristics deformation of weakly constrained riser in flow direction with different incident amplitudes of internal solitary waves

對于深海礦產輸送管道,輸送礦物濃度的變化會導致內部流體密度的變化,因此本文討論了內部流體不同密度對管道響應特征的影響.不同內部流體濃度條件下管道變形特征如圖12 所示.從圖中發現,管道內部流體密度對弱約束管道的流向變形影響較小.

圖12 內部流體密度對弱約束管道位移的影響Fig.12 Influence of inner fluid densities on displacement of weakly constrained riser

4 結論

本文基于勢流理論采用多域邊界單元法建立了內孤立波流場的計算模型,內孤立波波面采用全非線性動力學和運動學邊界條件,準確獲得了內孤立波流場的速度分布特征.上下層的流體質點速度方向相反,對于下凹型內孤立波,上層流體速度與波的傳播方向相同,下層速度與之相反,且上層流體質點速度明顯大于下層.

根據計算獲得的內孤立波流場信息,采用莫里森方程確定了內孤立波對海洋立管的載荷輸入,計算了頂張力立管和弱約束立管在內孤立波作用下的動力響應特征.管道的順流向位移和管道應力隨內孤立波波幅的增大而顯著,同時頂張力通過改變幾何剛度陣而對管道的動力響應產生明顯影響.由于上層流體速度明顯大于下層,且在所研究問題中拖曳力遠大于慣性力,因此管道順流向的最大位移發生在上層區域.同時計算了不同波幅和不同內部流體密度條件下弱約束管道的動力響應,發現內孤立波波幅對弱約束管道的動力響應影響顯著,而管道內部流體密度的影響則相對較小.

本文中內孤立波對立管的作用載荷根據流場信息采用莫里森方程給定,沒有考慮管道的渦激振動響應特征.為了更加全面分析管道在內孤立波作用下的動力響應規律,需要在今后的研究中進一步改進,為實際工程中的海洋立管設計和安全性評估提供更加充分的依據.