波流作用下柔性立管載荷響應數值分析

吳梓鑫， 陳昆鵬(上海船舶運輸科學研究所 航運技術與安全國家重點實驗室， 上海 200135)

波流作用下柔性立管載荷響應數值分析

吳梓鑫， 陳昆鵬
(上海船舶運輸科學研究所 航運技術與安全國家重點實驗室， 上海 200135)

結合商用軟件Fluent 14.5自定義函數(UDF)功能，采用速度邊界造波法和阻尼消波技術，通過流體體積函數(Volume of Fluid, VOF)追蹤自由液面，建立波浪、波流數值水池。采用雙向耦合分離法思想，借助Workbench 14.5流固耦合模塊System Coupling將流域和結構域計算數據交互傳遞，在數值波浪、波流水池中對長徑比為8的柔性立管進行雙向流固耦合模擬，分析其載荷響應，得出波浪場中不同流速下受力點的變化情況及不同波流方向下振動平衡狀態的變化情況。

數值水池； 流固耦合； 柔性立管； 載荷響應

0 引 言

目前國內外眾多科技工作者已在波浪與結構物相互作用的試驗研究和數值分析上取得許多成果，諸多結論已投入工程應用。例如：TAYLOR等[1]在小流速假定下用零流速時的三維脈動源加上流速修正項代替移動脈動源，利用繞射理論模擬三維浮體與波流的相互作用；周正全等[2]給出波流聯合作用下繞射問題的時域解理論模型，證明利用入射波勢和輔助輻射勢計算三維結構物波浪繞射力的可行性；SARPKAYA等[3-4]分析圓柱體阻力系數在波流場中的變化；HAN等[5]運用Hamiliton原理和Kirchhoff假設研究上端自由而下端支撐在海底條件下立管在波流作用下的軸向與橫向振動特性。

在數值研究方面：ATADAN[6]考慮基于非線性彈性理論的剪切效應，通過理論求解和數值模擬的方法確定海洋立管在波流作用下的動力響應；SHAH等[7]將海洋立管下端支撐而上端設置為簡支，研究其在波流聯合作用下的共振問題，并分析了流引起的主共振和波流聯合作用引起的組合共振特性；李軍強等[8]基于ANSYS和ABAQUS/Aqua等有限元分析軟件對波流聯合作用下深水立管的非線性動力響應進行諸多研究。

隨著商用軟件功能日趨完善，數值仿真已成為研究海洋問題的新方法，這里借助軟件Workbench 14.5模擬立管在數值水池中的運動響應，并進行相關分析。

1 數學模型

1.1 流體控制方程

假設流體為不可壓縮黏性流體，依據流體力學理論，其運動控制方程如下。

(1) 連續性方程為

(1)

式(1)中:ρ為流體密度;v為流場速度矢量。由于流體不可壓縮，連續性方程簡化為

(2)

(2) 對于黏性流體運動，當黏度為常數時，有動量守恒(Navier-Stokes，N-S)方程

(3)

ρ為常數時，流體不可壓縮，N-S方程簡化為

▽2v

(4)

1.2 結構控制方程

假定結構為線彈性，相對于其平衡位置作剛體運動和變形。結構經有限元離散后的動力學方程為

(5)

1.3 邊界條件

數值模擬中，流體計算主要需滿足自由液面條件和壁面條件。對柔性立管的計算還需滿足流固耦合界面條件。

1.3.1 自由液面條件

數值水池中流體上方為空氣，在自由液面處必須同時滿足運動學條件和動力學條件。對于運動學條件，針對不同的液面追蹤方法，有不同的數學表達方式;對于動力學條件，則須滿足表面應力條件。假設不計自由表面張力，表面應力條件可表示為

(6)

(7)

式(6)和式(7)中：un為垂直于自由表面的法向速度(流域中外法向為正)；ut為切向速度；p0為表面大氣壓強。

1.3.2 壁面條件

這里數值模擬中采用不可滑移壁面條件

v=vb

(8)

式(8)中：vb為相對于連體坐標系的運動速度。

1.3.3 流固耦合界面條件

流固耦合界面上邊界條件包括運動學條件和動力學條件。

(1) 運動學條件：流固交界面上法相速度保持連續。

vfn=vf·nf=vs·nf=-vs·ns=vsn

(9)

式(9)中：nf為流體邊界外法線向量；ns為固體邊界外法線向量。

(2) 動力學條件：流固交界面上法向應力保持連續。

σijnsj=τijnfj=-τijnsj

(10)

式(10)中：τij為流體應力張量的分量。

2 柔性立管在波流作用下載荷與響應分析

數值模擬借用軟件平臺Workbench 14.5，流體域通過Fluent計算，結構域選用Transient Structural計算，數據傳遞通過System Coupling模塊實現，后處理由CFD-Post完成。

2.1 模型建立及參數設置

2.1.1 流體域建模與設置

所模擬的余弦波波高為1 m，水池水深為5 m，周期為3.737 s，具體參數見表1。數值水池工作區60 m，消波區20 m，水深5 m，水面以上3 m為空氣，立管直徑1 m，距水池前段30 m，距左右邊界均為5 m。水池模型示意見圖1。

表1 余弦波參數

圖1 數值水池模型示意

采用結構網格劃分模型，在波長方向單個網格尺寸為1%波長，自由液面上下1個波高內網格高度為5%波幅，遠離自由液面處網格漸稀。以立管為中心，在xy平面立管周圍10 m×10 m內網格加密。

左端入口為速度入口，右邊界設置為壓力出口，上邊界設置為壓力入口，底邊界、水池兩壁面及立管邊界為無滑移固壁條件。造波、消波、壓力出口及自由液面捕捉通過Fluent軟件自定義函數(UDF)自編程序實現。

采用Fluent中三維非定常分離隱式求解器求解，湍流模型選用RNGk-ε，壓力方程選用加權體積力格式，壓力速度耦合方式采用PISO算法，體積分數方法為幾何重構，動網格更新方法選用擴散光順。流體域中與立管接觸面設置為流固耦合面。時間步長為0.005 s。

2.1.2 結構域建模與設置

立管直徑為1 m，壁厚為0.05 m，密度為7 850 kg/m3，具體參數見表2。立管與流域接觸面為流固耦合面，根據文獻[9]將底端設置為簡支，頂端限制其水平位移。此外，考慮到文獻[10]中涉及渦激振動和對預緊力進行計算分析，設置立管頂部預緊力為20 000 N。時間步長為0.005 s。立管表面網格與其接觸流體面網格并不需要完全對應，經過試算，單位網格尺寸設置為0.1 m。

2.1.3 流固耦合模塊設置

流體域和固體域分別設置完成后，通過System Coupling模塊將兩者同步連接，定義上述2部分流固耦合面為數據傳輸面，設定計算順序為先計算流體域再計算固體域，設置時間步長為0.005 s，計算時間為30 s。

2.2 數值模擬結果分析

為研究不同方向波流聯合作用下柔性立管的載荷響應，模擬表3中的工況。

表2 立管模型參數

表3 波流速度

為兼顧準確性和工程需要，采用文獻[11]中提出的波流水平速度疊加法進行波流聯合作用的數值模擬，以探求波浪參數在海流存在條件下的變化情況，其速度關系為

uwc(x,z,t)=uw(x,z,t)+uc

(11)

式(11)中：uwc(x,z,t)為波流聯合場水平速度；uw(x,z,t)為波浪場水平速度；uc為均勻流場水平速度。

2.2.1 立管受力分析

數值水池中波浪自左向右傳播，在立管迎浪面最前端z=3 m至z=-5 m每隔1 m設置1個點，共設置9個點監測其應力，壓力<0，拉力>0，立管在各水平高度下的壓力變化時程曲線見圖2。

a) Case 1

b) Case 2

c) Case 3

d) Case 4

在單獨波浪場中，立管最大受力位于立管底端，整體沿水平高度向上遞減，立管頂部承受拉力，各點受力周期相同。在波流聯合作用下，波流同向時Case 2中立管受力曲線分布與單獨波浪場中相似，最大受力位于立管底端，但在流的作用下立管受力上移。Case 3中立管受力分布發生顯著變化，立管受力向中間段趨近并以該水平高度向兩端逐漸減小，底端處由持續受壓狀態轉為壓力和拉力交替循環。當波流相向時，Case 4中立管的受力狀態與Case 2相近。

柔性立管在數值水池中的水平作用力系數時程曲線見圖3。波流同向時，在同一周期中由于流的作用，水平作用力系數較單獨波浪中正向幅值增大、負向幅值減小，且由于拖曳力正比于流速平方，Case 4中當流速增加時幅值增減變化更為明顯,當波流相向時水平作用力系數正向幅值正值比單獨波浪場中小，負向幅值比單獨波浪場大。

2.2.2 迎浪向振動分析

監測點在迎浪向的振動時程曲線見圖4，隨著波浪的傳遞作用，在約5 s時立管開始有較為明顯的振動，振動周期與波浪周期相近。波流同向時，由于流的作用加上立管兩端的約束和靜水壓力影響，在波浪傳至立管處前時監測點在迎浪向已有一定的位移，且在立管中間部分(即水深z=-1 m和z=-2 m區間內)位移最大，當波浪傳至立管處時監測點開始出現周期性振動，周期仍與立管單獨波浪中的振動周期基本相同。

圖3 柔性立管水平作用力系數曲線

Case 2中運動趨勢與單獨波浪中相似；Case 3中監測點振動曲線呈現出明顯的“峰尖谷坦”特性。波流相向時，由于反方向流的作用加上立管兩端約束條件的限制，立管迎浪面的初始位移有正有負，約為壓縮的“s”形。在波流聯合作用下，向迎浪向運動時速度較慢，向流向運動時速度較快，運動周期同樣為3.73 s。

各算例中柔性立管監測點在迎浪向的振幅匯總見圖5。雖然加載流的大小和方向存在差異，但振幅并沒有顯著的區別。穩定狀態時立管在z=0 m和z=-1 m處振幅最大，向上/下兩端逐漸減小，最大振幅區域均位于自由液面以下2 m內。

以流速和波速的比值作為橫坐標，選取算例中最大迎浪向振幅值(即各算例中z=-1 m處的振幅值)作為縱坐標，繪制成最大振幅趨勢曲線(見圖6)。由此可推斷出，在該模擬范圍內，當波流相向時，隨著流速減小，立管迎浪向振幅減小；當波流同向時，隨著流速增大，立管振幅增大，當流速與波速的比值約為25%時，振幅達到最大，此后流速繼續增大，立管振幅減小。

a) Case 1

b) Case 2

c) Case 3

d) Case 4

圖5 監測點迎浪向振幅匯總

圖6 立管迎浪向最大振幅隨波流比值走勢

立管在數值水池中振動時其平衡位置相對于靜止位置發生偏移，各算例立管在迎浪向振動的平衡位置偏移原靜止位置(見圖7)。在單獨波浪場中，立管振動平衡位置較靜止位置在迎浪向有一定的位移，最大振動平衡位置區域位于自由液面以下約4 m處。加載流之后，波流同向時，最大平衡位置上移至自由液面以下1～2 m區域內，當流速為波速的20%時最大偏移值約為單獨波浪場最大偏移值的8倍，流速加大到40%波速時最大偏移值約為單獨波浪場中23倍；波流相向時，各水平位置的振動平衡位置的變化與波流同向時差異較大。由于流的作用對立管振動平衡位置的影響非常顯著，因此在波流相向時，立管背浪面自由液面以下在流等外載荷作用下其振動平衡位置向迎浪向負方向偏移；而立管結構具有各向同性的屬性使得迎浪面自由液面以下其平衡位置向迎浪向偏移，立管頂端水平位移受限，因此其自由液面以上部分振動平衡位置在原靜止位置左側，Case 4中立管迎浪面和背浪面立管振動平衡位置見圖8。

圖7 立管迎浪向振動平衡位置

圖8 Case 4立管振動平衡位置

2.2.3 橫浪向振動分析

數值水池中立管監測點在橫浪向的振動時程曲線見圖9。相比立管迎浪向振動，橫浪向振動隨機性更強，振動周期差異較大，但振幅最大區域依然位于自由液面以下2 m區域內。分析數據可得，當只有波浪作用時，橫浪向振幅不到迎浪向振幅的0.1%。當波流共同作用時，橫浪向振動不可忽略，波流同向中流速為波速20%時橫浪向振幅為迎浪向振幅的17.18%，當流速增至波速的40%時橫浪向振幅為迎浪向振幅的24.2%；而波流相向時，橫浪向振幅為迎浪向振幅的1.11%。

a) Case 1

b) Case 2

c) Case 3

d) Case 4

2.2.4 軸向振動分析

監測點在立管軸向的振動時程曲線見圖10。與水平面內振動不同，由于監測點在水平面內具有運動位移，而立管底端簡支約束，上端只能進行軸向運動，因此各點在軸向的運動方向不再統一，當靠近底端部分位移>0時頂端位移<0，反之亦然。幅值方面，在立管頂部軸向運動幅值最大，而在自由液面以下2 m區域內軸向運動幅值最小。在單獨波浪場和正負向流速為波速20%的波流數值水池中，各點軸向運動時程曲線約為余弦函數狀，而當流速增至波速的40%時，時程曲線在單個周期內存在2個正向幅值和2個負向幅值，體現出流速較大時波流聯合作用下立管運動的復雜性。此外，數據分析表明，立管軸向運動的周期與迎浪向的運動周期基本吻合。

a) Case 1

b) Case 2

c) Case 3

d) Case 4

圖11 監測點軸向振幅

柔性立管在軸向運動振幅統計見圖11。各算例中立管各水平位置在軸向的振幅變化趨勢一致，均在立管中間位置，即迎浪向振幅最大位置z=-1 m處振幅最小，向上和向下振幅都逐漸增大。由于底端簡支約束，振幅為0，且振幅由z=-1 m處向下為0的中間部分存在1個幅值轉捩點，由圖11可知該轉捩點在z=-3 m附近。

3 結 語

基于黏流理論，以商業軟件Fluent 14.5為平臺，通過程序接口加載自編程序進行開發，采用速度邊界條件造波法建立效果良好的波流數值水池。借助軟件平臺Workbench 14.5，通過流固耦合模塊System Coupling連接流體模塊Fluent和結構模塊Transient Structural，數值模擬長徑比為8的柔性立管在不同波流參數數值水池中的運動，從立管受力、迎浪向、橫浪向和軸向振動等方面對柔性立管在波流聯合作用下的載荷響應進行分析和研究。結果表明:在波浪場中，立管底端受力最大，沿軸向向上逐漸減?。划敿虞d流速較小時，立管受力分布變化不大；當流速較大時，立管應力集中點上移。在模擬范圍內，波流同向時，立管振動平衡位置向迎浪向偏移，迎浪向振幅隨流速先增大，達到某一極值后減??；波流相向時，振幅大于單獨波浪場中的振幅。此外，單獨波浪場中，立管橫浪向振幅很小，而在波流聯合作用下橫浪向振幅顯著增大，不可忽略。

[1] TAYLOR E R, HU C S, NIELSEN F G. Mean Drift Forces on Slowly Advancing Vertical Cylinders in Long Waves[J]. Applied Ocean Research,1990,12(3):141-152.

[2] 周正全，張亮，戴遺山. 波流聯合作用下的物體上的波浪繞射力[J]. 水動力學研究與進展 (A 輯)，1992，7(4)：476-483.

[3] SARPKAYA T, BAKMIS C, STORM M A. Hydrodynamic Forces from Combined Wave and Current Flow on Smooth and Rough Circular Cylinders at High Reynolds Number[C].Proceeding of 16thOTC, Houston,1984.

[4] SARPKAYA T, STORM M. In-Line Force on a Cylinder Translating in Oscillatory Flow[J]. Applied Ocean Research, 1985,7(4):188-196.

[5] HAN S M, BENAROYA H. Non-Linear Coupled Transverse and Axial Vibration of a Compliant Structure, Part 1: Formulation and Free Vibration[J]. Journal of Sound and Vibration, 2000, 237(5): 837-873.

[6] ATADAN A S, CALISAL S M, MODI V J, et al. Analytical and Numerical Analysis of the Dynamics of a Marine Riser Connected to a Floating Platform[J]. Ocean Engineering, 1997, 24(2): 111-131.

[7] SHAH K, FERZIGER J H. A Fluid Mechanician's View of Wind Engineering: Large Eddy Simulation of Flow Past a Cubic Obstacle[J].Wind Engineering, 1997, 67(4): 211-224.

[8] 李軍強，劉宏昭，何欽象,等. 波浪力作用下海洋鉆井隔水管隨機振動研究[J]. 機械科學與技術，2004，23(1): 7-10.

[9] 羅冬冬. 剪切流下長徑比對柔性立管渦激振動的影響[D].鎮江：江蘇科技大學，2014.

[10] 王成官. 海洋隔水管渦激振動特性的三維數值模擬研究[D].上海：上海交通大學，2011.

[11] 李玉成. 波浪與水流共同作用下的流速場[J]. 海洋工程，1983(4):12-23.

Test of Loads on Flexible Riser and its Response with Numerical Wave-Current Tank

WUZixin,CHENKunpeng
(State Key Laboratory of Navigation and Safety Technology, Shanghai Ship & Shipping Research Institute, Shanghai 200135,China)

Numerical wave/wave-current tanks with high effectiveness in wave-making and wave-absorbing functions are established on the platform of Fluent 14.5 with UDF programs by using velocity boundary and damping wave technique as well as VOF method. With the two-way coupling separation concept, transferring the data of the fluid domain and structure domain in the module System Coupling of the Workbench 14.5, the motion response of a flexible riser with the aspect ratio of 8 is simulated in the numerical tank, and the loads and response on the riser are analyzed.

numerical tank; Fluid-Structure Interaction; flexible riser; loads and response

2016-08-30

吳梓鑫(1989—)，男，江蘇南通人，研究實習員，碩士，主要從事船舶水動力研究。

1674-5949(2016)04-0014-08

P756.2

A