波流作用下柔性立管載荷響應(yīng)數(shù)值分析

吳梓鑫， 陳昆鵬(上海船舶運(yùn)輸科學(xué)研究所 航運(yùn)技術(shù)與安全國家重點(diǎn)實(shí)驗(yàn)室， 上海 200135)

波流作用下柔性立管載荷響應(yīng)數(shù)值分析

吳梓鑫， 陳昆鵬
(上海船舶運(yùn)輸科學(xué)研究所 航運(yùn)技術(shù)與安全國家重點(diǎn)實(shí)驗(yàn)室， 上海 200135)

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

數(shù)值水池； 流固耦合； 柔性立管； 載荷響應(yīng)

0 引 言

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

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

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

1 數(shù)學(xué)模型

1.1 流體控制方程

假設(shè)流體為不可壓縮黏性流體，依據(jù)流體力學(xué)理論，其運(yùn)動(dòng)控制方程如下。

(1) 連續(xù)性方程為

(1)

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

(2)

(2) 對于黏性流體運(yùn)動(dòng)，當(dāng)黏度為常數(shù)時(shí)，有動(dòng)量守恒(Navier-Stokes，N-S)方程

(3)

ρ為常數(shù)時(shí)，流體不可壓縮，N-S方程簡化為

▽2v

(4)

1.2 結(jié)構(gòu)控制方程

假定結(jié)構(gòu)為線彈性，相對于其平衡位置作剛體運(yùn)動(dòng)和變形。結(jié)構(gòu)經(jīng)有限元離散后的動(dòng)力學(xué)方程為

(5)

1.3 邊界條件

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

1.3.1 自由液面條件

數(shù)值水池中流體上方為空氣，在自由液面處必須同時(shí)滿足運(yùn)動(dòng)學(xué)條件和動(dòng)力學(xué)條件。對于運(yùn)動(dòng)學(xué)條件，針對不同的液面追蹤方法，有不同的數(shù)學(xué)表達(dá)方式;對于動(dòng)力學(xué)條件，則須滿足表面應(yīng)力條件。假設(shè)不計(jì)自由表面張力，表面應(yīng)力條件可表示為

(6)

(7)

式(6)和式(7)中：un為垂直于自由表面的法向速度(流域中外法向?yàn)檎?；ut為切向速度；p0為表面大氣壓強(qiáng)。

1.3.2 壁面條件

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

v=vb

(8)

式(8)中：vb為相對于連體坐標(biāo)系的運(yùn)動(dòng)速度。

1.3.3 流固耦合界面條件

流固耦合界面上邊界條件包括運(yùn)動(dòng)學(xué)條件和動(dòng)力學(xué)條件。

(1) 運(yùn)動(dòng)學(xué)條件：流固交界面上法相速度保持連續(xù)。

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

(9)

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

(2) 動(dòng)力學(xué)條件：流固交界面上法向應(yīng)力保持連續(xù)。

σijnsj=τijnfj=-τijnsj

(10)

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

2 柔性立管在波流作用下載荷與響應(yīng)分析

數(shù)值模擬借用軟件平臺(tái)Workbench 14.5，流體域通過Fluent計(jì)算，結(jié)構(gòu)域選用Transient Structural計(jì)算，數(shù)據(jù)傳遞通過System Coupling模塊實(shí)現(xiàn)，后處理由CFD-Post完成。

2.1 模型建立及參數(shù)設(shè)置

2.1.1 流體域建模與設(shè)置

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

表1 余弦波參數(shù)

圖1 數(shù)值水池模型示意

采用結(jié)構(gòu)網(wǎng)格劃分模型，在波長方向單個(gè)網(wǎng)格尺寸為1%波長，自由液面上下1個(gè)波高內(nèi)網(wǎng)格高度為5%波幅，遠(yuǎn)離自由液面處網(wǎng)格漸稀。以立管為中心，在xy平面立管周圍10 m×10 m內(nèi)網(wǎng)格加密。

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

采用Fluent中三維非定常分離隱式求解器求解，湍流模型選用RNGk-ε，壓力方程選用加權(quán)體積力格式，壓力速度耦合方式采用PISO算法，體積分?jǐn)?shù)方法為幾何重構(gòu)，動(dòng)網(wǎng)格更新方法選用擴(kuò)散光順。流體域中與立管接觸面設(shè)置為流固耦合面。時(shí)間步長為0.005 s。

2.1.2 結(jié)構(gòu)域建模與設(shè)置

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

2.1.3 流固耦合模塊設(shè)置

流體域和固體域分別設(shè)置完成后，通過System Coupling模塊將兩者同步連接，定義上述2部分流固耦合面為數(shù)據(jù)傳輸面，設(shè)定計(jì)算順序?yàn)橄扔?jì)算流體域再計(jì)算固體域，設(shè)置時(shí)間步長為0.005 s，計(jì)算時(shí)間為30 s。

2.2 數(shù)值模擬結(jié)果分析

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

表2 立管模型參數(shù)

表3 波流速度

為兼顧準(zhǔn)確性和工程需要，采用文獻(xiàn)[11]中提出的波流水平速度疊加法進(jìn)行波流聯(lián)合作用的數(shù)值模擬，以探求波浪參數(shù)在海流存在條件下的變化情況，其速度關(guān)系為

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

(11)

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

2.2.1 立管受力分析

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

a) Case 1

b) Case 2

c) Case 3

d) Case 4

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

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

2.2.2 迎浪向振動(dòng)分析

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

圖3 柔性立管水平作用力系數(shù)曲線

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

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

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

a) Case 1

b) Case 2

c) Case 3

d) Case 4

圖5 監(jiān)測點(diǎn)迎浪向振幅匯總

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

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

圖7 立管迎浪向振動(dòng)平衡位置

圖8 Case 4立管振動(dòng)平衡位置

2.2.3 橫浪向振動(dòng)分析

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

a) Case 1

b) Case 2

c) Case 3

d) Case 4

2.2.4 軸向振動(dòng)分析

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

a) Case 1

b) Case 2

c) Case 3

d) Case 4

圖11 監(jiān)測點(diǎn)軸向振幅

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

3 結(jié) 語

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

[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] 周正全，張亮，戴遺山. 波流聯(lián)合作用下的物體上的波浪繞射力[J]. 水動(dòng)力學(xué)研究與進(jìn)展 (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] 李軍強(qiáng)，劉宏昭，何欽象,等. 波浪力作用下海洋鉆井隔水管隨機(jī)振動(dòng)研究[J]. 機(jī)械科學(xué)與技術(shù)，2004，23(1): 7-10.

[9] 羅冬冬. 剪切流下長徑比對柔性立管渦激振動(dòng)的影響[D].鎮(zhèn)江：江蘇科技大學(xué)，2014.

[10] 王成官. 海洋隔水管渦激振動(dòng)特性的三維數(shù)值模擬研究[D].上海：上海交通大學(xué)，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—)，男，江蘇南通人，研究實(shí)習(xí)員，碩士，主要從事船舶水動(dòng)力研究。

1674-5949(2016)04-0014-08

P756.2

A