基于CFD方法的FPSO的波浪荷載和運動響應研究?

左伊人,周 琦,楊風艷,張 敏??

(1. 中國海洋大學工程學院,山東省海洋工程重點實驗室,山東 青島 266100;2. 浙江華東建設工程有限公司,浙江 杭州 310014; 3. 海洋石油工程(青島)有限公司,山東 青島 266520)

浮式生產儲油船(Floating production storage and offloading,FPSO)是目前海上油氣田開發生產的主流設施之一。隨著海洋工程向深水推進,工程所處的海洋環境條件愈加惡劣,體現在波浪出現強烈的非線性,而波浪對其船體結構施加的荷載和造成的運動響應是工程上需要重點關注的問題。為保證FPSO在服役時的安全性和工作效率,對于其波浪荷載的計算以及運動響應的預報必須要求準確。

對于FPSO所受的波浪荷載,早期多是使用勢流理論進行計算,Newman[1]提出的遠場法以及Pinkster等[2]提出的近場法均得到廣泛的應用,但是勢流理論較難考慮流體黏性效應和強非線性效應,不能很好地反映實際工程面臨的真實問題。Inoue等[3]使用近場法和遠場法計算了液化天然氣浮式生產儲油船 (Liquefied natural gas FPSO,LNG-FPSO)系統的漂移力,同模型試驗對比發現,不考慮黏性會導致漂移力的預測誤差。國內的學者王科等[4]指出計算FPSO的縱搖和橫搖需要考慮黏性阻尼。隨著計算機技術和數值方法的快速發展,計算流體力學(Computational fluid dynamics,CFD)方法以其考慮流體黏性和非線性影響的優勢開始受到廣泛關注。外國學者Huijsmans[5]首次使用CFD方法計算了深水規則波作用下FPSO的二階波浪力的傳遞函數。Rosetti等[6]預測了不同波陡的規則波作用下二維FPSO的波浪力,發現CFD方法可以很好地捕捉到力的最大值。Hu等[7]使用CFD求解器計算了極端波作用下FPSO的力和運動,且通過同勢流結果和實驗數據的對比發現,縱搖的黏性效應明顯。王和靜等[8]使用CFD方法計算了FPSO的平均漂移力和運動響應,發現較之勢流,CFD方法會高估平均漂移力。李奇等[9]應用CFD求解器naoe-FOAM-SJTU對單點系泊的FPSO水動力性能進行數值模擬,驗證了其準確性,體現出CFD進行強非線性現象仿真的優勢。但是,強非線性波浪力的作用機理目前還不明晰,未提出有效的強非線性波浪荷載的預報模型[10]。

本文應用CFD軟件STAR-CCM+[11]建立三維數值波浪水池,使用Richardson外推法分析了計算收斂性和數值不確定度。在數值水池中模擬了規則波作用下FPSO的垂蕩與縱搖,計算迎浪波浪力,分析流場瞬時自由液面和渦量場。通過同勢流結果對比,探究了黏性效應對波浪力和運動響應的影響,探究了非線性波浪力的作用機理。本文的數值模擬使用實尺度FPSO,避免船體尺度效應,充分考慮了波浪與FPSO相互作用過程中的黏性效應,對計算及修正浮體受到的二階波浪力和運動響應具有一定的理論指導意義。

1 數值分析方法

本文以雷諾平均方程(RANS equation)為控制方程求解瞬態的黏性繞流場,采用有限體積法(Finite volume method,FVM)離散控制方程。選擇歐拉多相流創建水和空氣兩相流體,流域為不可壓縮流體,采用流體體積法(Volume of fluid,VOF)捕捉自由液面,用SSTk-ω湍流模型封閉方程,使用壓力耦合方程組的半隱式方法(Semi-implicit method for pressure linked equations,SIMPLE)耦合求解速度和壓力,并考慮重力的影響。對于浮體的運動,采用重疊網格技術和動態流體固體相互作用模型DFBI來模擬。

連續方程:

(1)

RANS方程:

(2)

本文采用SSTk-ω湍流模型包含的k方程和ω方程分別為

(3)

(4)

式中:k為湍動能;ω為耗散率;Γk和Γω是有效擴散項;Gk和Gω是生成項;Yk和Yω是湍動耗散項;Dω是交叉擴散項;Sk和Sω是自定義源項。

自由液面用VOF法處理,體積占比aq滿足方程:

(5)

2 數值模擬

2.1 船型幾何參數

本研究采用的船型為FPSO[8],采用實船尺寸以避免船舶尺度效應。計算模型參數見表1。后文中的船長均指垂線間長Lpp。

表1 計算模型參數Table 1 Parameters of the FPSO

2.2 計算域設置

本文根據文獻[12]的建議對計算域進行設置:整個計算域長度約為船長的9.8倍,寬度約為船長的6.4倍;速度入口距離船首為船長的3.2倍,壓力出口距離船尾為船長的6.6倍,且整個計算域包含2倍波長的消波區域;整個場域的上部為空氣,下部為水,船下部邊界距船體為船長的2倍。FPSO在數值水池中的設置如圖1所示,因船體對稱,使用FPSO一半船體進行計算。原點為O,橫向(y軸)位置在中線面上,縱向(x軸)位置在船中,豎向(z軸)位置在距船舶基線18.44 m處。

圖1 數值水池示意圖Fig.1 Sketch of the numerical tank

2.3 造波和消波

采用邊界造波法原理,直接通過雷諾平均納維-斯托克斯(Reynolds-averaged Navier-Stokes equations,RANS)求解器分別控制速度入口和壓力出口參數來模擬造波。入口邊界、側面邊界設置為速度入口,出口邊界設置為壓力出口,3個邊界均設置強制消波[13],在消波區域強制讓波浪保持入射的狀態。船體表面為不可滑移固壁條件,近壁面采用壁面函數法處理。計算域頂部和底部均設置為速度入口,以更好地模擬波浪前進,減輕波浪衰減。整個流場造波效果見圖2,波面高程的原點在-11.45 m。

圖2 空場造波的自由液面高程圖Fig. 2 Free surface elevation of wave generation

2.4 時間步長及網格尺寸收斂性驗證

使用重疊網格,根據一個波高范圍內劃分10～20個網格,一個波長范圍內劃分80～120個網格的原則,網格劃分如圖3所示。

圖3 網格劃分的示意圖Fig.3 Sketch of the mesh

表2 網格尺寸方案Table 2 Test case of grid size

根據文獻[12]的建議,設置3種時間步長,如表3所示。

首先通過空場造波驗證數值波浪的可行性。3套網格下的造波結果與理論值的對比如圖4所示。由圖4可知,3套網格的數值結果同理論值匹配較好。

圖4 不同網格尺寸條件下的造波與理論值對比Fig.4 Comparison of the numerical and theoretical free surface elevation

(6)

(7)

(8)

基于Richardson外推法[14-15],可以定義收斂率R、精度P、外推值Sext、相對誤差估計值ea、相對誤差外推值eext和網格收斂指數(Grid convergence index,GCI)I,以對時間步長和網格尺寸進行收斂性研究,并分析數值不確定度:

(9)

(10)

(11)

(12)

(13)

(14)

式中:s1、s2、s3分別為在細密網格(或最小時間步長)、中等網格(或中等時間步長)和粗糙網格(或最大時間步長)計算得到的數值結果(運動響應、力);Sext,32表示基于s3、s2的外推值;ea,32表示基于s3、s2的相對誤差估計值;eext,32表示基于s3、s2的相對誤差外推值。

根據收斂率R的大小,可以分成4種收斂情況:單調收斂(01);振蕩發散(R<-1)。其中,單調收斂情況使用網格收斂指數表征不確定度,振蕩收斂情況使用數值最大值和最小值的差值的一半來表征不確定度。表4和5分別給出了時間步長收斂性和網格收斂性的基于s2的收斂性驗證結果。

表4 時間步長收斂性Table 4 Time-step convergence

根據表4可以發現,采用0.02 s的時間步長已經滿足收斂。采用最小時間步長0.015 s得到的計算結果數值不確定度很低,垂蕩的不確定度小于5%,縱搖和迎浪平均波浪力系數的不確定度更是小于1%。根據表5可以看出使用粗糙網格已經滿足收斂。使用中等網格的垂蕩和迎浪平均波浪力系數的不確定度均小于1%,且縱搖不確定度遠小于1%,說明網格收斂性很好。在滿足時間步長收斂和網格收斂的精度條件下,為了保證計算的物理時間效率,后續選擇0.015 s作為時間步長。

表5 網格收斂性Table 5 Grid convergence

2.5 數值模型驗證

使用本文的數值水池,選擇和文獻[8]同樣的波浪條件進行計算,時歷曲線對比的結果如圖5所示。結果吻合良好,證明本文的數值模型具有可行性。

圖5 數值計算結果與參考結果的對比圖Fig.5 Comparison between numerical results and reference results

2.6 計算工況

為更加真實地模擬波浪,計算中使用斯托克斯五階波 (Stokes fifth order wave)和0°浪向角(迎浪方向)。根據文獻[16]中的實際海域環境條件,分別選擇5組同一周期不同波高的工況和5組不同周期同一波高的工況(見表6)。由于不考慮系泊,計算過程僅開放縱搖和垂蕩兩個自由度。

表6 計算工況Table 6 Calculation conditions

3 計算結果及分析

計算所有工況的垂蕩響應、縱搖響應和迎浪波浪力。選取FPSO運動較為穩定的10個周期。

3.1 波浪荷載

圖6展示了在一個波浪周期內不同波高和不同波長工況下的迎浪波浪力時歷曲線。由圖6(a)可知,一個波浪周期內,每種波高的工況下均存在2個極大值和2個極小值,且2個極大值的幅值較為接近,2個極小值的幅值也很接近。由圖6(b)可知,在一個波浪周期內,僅在波長不大于船長的工況下可明顯看到波浪力存在第二個極值。

圖6 一個周期內迎浪波浪力的時歷曲線Fig.6 Time history of longitudinal wave forces in one wave period

對所有工況的迎浪波浪力時歷曲線做頻譜分析。由圖7(a)可知,在不同波高的工況下,0.18 Hz(二倍波頻)附近對應的力(二倍波頻力)均非常大,0 Hz(0倍波頻)附近對應的力(平均波浪力)也較為明顯,同一波浪頻率對應的力數值接近。0.27 Hz(三倍波頻)附近、0.36 Hz(四倍波頻)附近、0.45 Hz(五倍波頻)附近均有明顯數值,說明迎浪波浪力的非線性強烈。由圖7(b)可知,波長大于船長的2個工況下的二倍波頻對應的力非常小,其他3個工況下的二倍波頻對應的力均為波頻對應的力的一半以上。綜合圖6和7分析得知,波浪力時歷曲線中,一個周期內存在2個極大值和2個極小值,且均同二倍波頻有關,二倍波頻對應的力越大,第二個極大值和第二個極小值越明顯。

圖7 迎浪波浪力的頻譜圖Fig.7 Frequency spectrum of the longitudinal wave forces

圖8和9分別展示了平均波浪力和二倍波頻力隨波浪參數的變化趨勢,其中的勢流結果由AQWA(Advanced quantitative wave analysis)軟件[17]計算得到。從圖8可以看出平均波浪力和二倍波頻力隨波高增大呈指數增大。隨著波高增大,迎浪平均波浪力的CFD仿真結果同其勢流結果的誤差逐漸增大,而倍波頻力的CFD仿真結果同其勢流結果的差異不明顯。說明波高增大時,黏性效應對平均波浪力的影響明顯。從圖 9 可以看出,隨著波長L變化,波浪力變化較為復雜。隨波長增大,迎浪平均波浪力而減小,當波長L< 169.5 m 時變化速率較慢,且此時平均波浪力的 CFD結果同其勢流結果有誤差較大。隨波長增大,二倍波頻力先增大再減小;當波長接近船長時,二倍波頻力達到極大值;當波長L<169.5 m 時,二倍波頻力的 CFD 結果也同其勢流結果有較大誤差。說明當波長較短(L< 169.5 m)時,黏性效應對二階波浪力存在明顯影響。

圖8 二階波浪力隨波高的變化圖Fig.8 Variation of second-order wave forces with wave height

圖9 二階波浪力隨波長的變化圖Fig.9 Variation of second-order wave forces with wave length

3.2 FPSO的運動響應

圖10展示了波高變化的工況里一個波浪周期內FPSO的垂蕩和縱搖的時歷曲線,時歷曲線均呈現余弦函數的特征,且波高越大,運動幅值越大。

圖10 波高不同工況中一個波浪周期內垂蕩(a)和縱搖(b)的時歷曲線Fig.10 Time history of heave (a) and pitch (b) in one time period of different wave height

圖11為響應幅值算子(Response amplitude operator,RAO)隨波高的變化圖,其中,勢流結果不隨波高變化。圖12為RAO隨波長的變化圖,RAO變化復雜,短波時的變化速率較快,長波時的變化速率較慢。可以看出垂蕩的CFD結果與其勢流結果吻合良好。而在波高較大的情況下,縱搖的CFD結果同其勢流結果存在差異。

圖11 RAO隨波高的變化圖Fig.11 Variation of RAO with wave height

圖12 RAO隨波長的變化圖Fig.12 Variation of RAO with wave length

為定量分析黏性對縱搖的影響,本文中比較了縱搖的CFD結果和勢流結果的相對誤差(見表7)。小波陡情況下,CFD結果同勢流結果的相對誤差小于5%。當波陡H·L-1>1/15時,相對誤差達到了10%以上。說明波陡較大(H·L-1>1/20)時,黏性效應對于縱搖的影響明顯。

表7 縱搖運動的CFD結果和勢流結果的相對誤差Table 7 Relative error of pitch motion between CFD results and potential flow results %

3.3 瞬時時刻流場分析

圖13以工況5為例,展示了一個周期內FPSO在流場中的運動以及自由液面的演變和壓力分布的時刻圖。圖13中選取的時刻分別對應受力曲線極值點的時刻,目的在于直觀地觀察波浪入射和FPSO的運動,進而分析二者的相互影響對FPSO受力的影響。

圖13 工況5的自由液面和壓力分布時刻圖Fig.13 Snap shots of wave pattern and pressure of condition 5

t=152.4 s時,FPSO剛運動過平衡位置,船身下沉并尾傾,此時波峰位于船首和船尾,船身壓力分布均勻。船尾柱下沉過程有明顯的壓力集中。t=155.4 s時,FPSO的垂蕩達到最大幅值,船尾柱下沉并且船首傾,此時船首遭遇波谷,使得波谷的峰值更小。波峰位于船身前中段,壓力集中在此區域。t=158.4 s時,波峰即將傳播到船首,使壓力集中于船首并導致其下沉,而船中后段正處于波峰的位置,此區域壓力數值較大。t=161.4 s時,波峰傳播到船首,船身同時上浮,使得船首壓力數值非常大,而船中和船尾部的壓力均較小,進而使船體首、尾間的壓力差非常大。

本文使用的FPSO船首為直立U形結構,而船尾為傾斜方形結構且有尾柱,二者結構的差異對流場的擾動完全不同。圖13能看出船首前端的波峰和波谷都有明顯增大,而船尾的流場擾動不如船首明顯。圖14展示了同一時刻船首、尾處的垂向截面的渦量與流線圖,箭頭為流線的方向。t=150 s時,FPSO正在首傾,船首處的流體流動方向不同于船尾的流體流動方向,船尾柱近壁面能看到渦量,且流線復雜。

圖14 同一時刻船首、尾垂向截面的渦量場和流線圖Fig.14 The vorticity and streamlines of vertical section of the bow and stern at the same time

同一時刻不同波陡條件下的水下某水平截面的渦量場(底色)和流線(黑色線條)如圖15所示。隨著波陡增大,船首附近的流線同FPSO近壁面的渦量場也存在差異。

綜合分析,本文作者認為CFD方法可以良好地捕捉到流場的特征。FPSO在波浪中不斷運動,濕表面積和船體壓力不斷變化,從而導致自由液面的復雜擾動,影響了受力的非線性。在大波陡情況下,船首波浪疊加現象明顯,自由液面復雜,船首流場同船尾流場有明顯差異,導致黏性效應對波浪力和縱搖運動的影響顯著。

4 結論

(1)考慮黏性的CFD方法較之勢流理論會高估平均波浪力。波高越大或者波長越短,黏性效應越強,從會導致CFD結果和勢流結果的差值差異化增大。

(2)規則波作用下FPSO受力呈現非線性的原因:入射波與FPSO運動相互影響導致自由液面的復雜變化;FPSO運動時濕表面積和壓力分布不斷變化。

(3)由于船首結構不同于船尾結構,導致流場擾動出現明顯差異,當波陡達到1/15,CFD結果與勢流結果的誤差達到10%。說明在非線性較強的情況下,黏性效應對縱搖運動影響明顯。