基于MA TLAB的海洋立管渦激振動(dòng)數(shù)值模擬系統(tǒng)研究

李洪春,郭海燕,李效民

(中國(guó)海洋大學(xué)工程學(xué)院,山東青島266100)

基于MA TLAB的海洋立管渦激振動(dòng)數(shù)值模擬系統(tǒng)研究

李洪春,郭海燕,李效民

(中國(guó)海洋大學(xué)工程學(xué)院,山東青島266100)

基于MATLAB開(kāi)發(fā)海洋立管渦激振動(dòng)數(shù)值模擬系統(tǒng)NSVIV 1.0,系統(tǒng)采用尾流振子模型模擬外流對(duì)立管結(jié)構(gòu)的作用,考慮內(nèi)流對(duì)立管結(jié)構(gòu)的影響,對(duì)立管的渦激振動(dòng)動(dòng)力響應(yīng)和疲勞壽命進(jìn)行預(yù)測(cè)分析。系統(tǒng)界面簡(jiǎn)潔清晰,使用方便,適用于頂張力立管,為進(jìn)一步集成功能齊全的海洋立管設(shè)計(jì)分析軟件打下基礎(chǔ)。

MATLAB;圖形用戶(hù)界面;海洋立管;渦激振動(dòng);數(shù)值模擬

立管是海面作業(yè)平臺(tái)與海底井口間的主要連接件,是海洋油氣開(kāi)發(fā)的重要設(shè)施。立管運(yùn)行期間受力復(fù)雜,當(dāng)波浪或海流流經(jīng)立管時(shí),在一定的流速條件下,可在立管兩側(cè)交替形成旋渦,當(dāng)旋渦脫落頻率與立管固有頻率接近時(shí),將使立管強(qiáng)烈振動(dòng),發(fā)生“鎖定”(lock-in)現(xiàn)象,加劇立管的疲勞破壞。因此,渦激振動(dòng)引起的疲勞損傷,是立管設(shè)計(jì)中必須考慮的重要因素。

立管渦激振動(dòng)的數(shù)值模擬研究,主要有2種方法[1]:一種是通過(guò)求解Navier-Strokes方程的計(jì)算流體動(dòng)力學(xué)法(CFD)。另一種是通過(guò)實(shí)驗(yàn)數(shù)據(jù)確定流體力參數(shù)的半經(jīng)驗(yàn)法,基于此種方法的幾種計(jì)算模型在海洋工程中得到廣泛的應(yīng)用,比較著名的有MIT的Vandiver教授開(kāi)發(fā)的SHEAR7[2]、MIT的Triantafyllou教授開(kāi)發(fā)的VIVA[3]、挪威的Larsen等開(kāi)發(fā)的VIVANA[4]等。

內(nèi)流對(duì)立管振動(dòng)的影響已引起愈來(lái)愈多學(xué)者的重視,郭海燕等[5-7]考慮管內(nèi)流體和管外海洋環(huán)境荷載的共同作用,分析立管的靜、動(dòng)力特性,并開(kāi)展了相關(guān)理論和實(shí)驗(yàn)研究。

目前國(guó)外的渦激振動(dòng)分析軟件未考慮立管內(nèi)部流體的作用,本文開(kāi)發(fā)的海洋立管渦激振動(dòng)數(shù)值模擬系統(tǒng)—NSVIV 1.0,考慮內(nèi)流作用,采用半經(jīng)驗(yàn)的尾流振子模型模擬流體力,用時(shí)程分析法,求解立管在任意剖面形式海流作用下的渦激振動(dòng)響應(yīng)和疲勞壽命,計(jì)算結(jié)果將更接近工程實(shí)際。系統(tǒng)基于常用數(shù)學(xué)軟件MATLAB開(kāi)發(fā),界面簡(jiǎn)潔清晰,使用方便。NSVIV 1.0適用于頂張力立管的分析,隨著研究的不斷深入,將升級(jí)為包括鋼懸鏈線立管等在內(nèi)、功能齊全的海洋立管渦激振動(dòng)分析軟件,最終有望用于我國(guó)海洋立管的工程設(shè)計(jì)。

1 計(jì)算原理

1.1 立管動(dòng)力學(xué)模型

立管受內(nèi)部流體和外部海洋環(huán)境荷載的共同作用,見(jiàn)圖1。假定管外海流沿X軸正向流動(dòng),管內(nèi)流體自上而下流動(dòng)。立管將在YOZ平面內(nèi)產(chǎn)生橫向渦激振動(dòng),在XOZ平面內(nèi)產(chǎn)生順流向振動(dòng),本系統(tǒng)主要研究立管橫向的動(dòng)力響應(yīng),取立管單元dz進(jìn)行分析,見(jiàn)圖2。

圖1 海洋立管振動(dòng)示意圖Fig.1 Vibration of marine riser

圖2 立管單元dz變形圖Fig.2 Deformation of riser element dz

根據(jù)功-能平衡原理,可得立管橫向渦激振動(dòng)方程[8]:

其中:m=mr+mi+ma,mr為單位長(zhǎng)度立管質(zhì)量,mi為單位長(zhǎng)度管內(nèi)流體質(zhì)量,ma為附加流體質(zhì)量;y為立管橫向位移;c為結(jié)構(gòu)阻尼系數(shù),c′為流體阻尼系數(shù);V為內(nèi)部流體流速;E為立管彈性模量,I為立管截面慣性矩;Te表示立管單元的有效軸向力;ρsea為海水密度,U為外流流速;D為立管直徑;CL0為升力系數(shù)幅值,一般取CL0=0.3;q為升力振子系數(shù),q=2CL/CL0。

1.2 流體力模型

本文采用Facchinetti et al.[9]推薦的尾流振子模型,該模型只考慮流體與結(jié)構(gòu)之間的線性耦合項(xiàng),其表達(dá)式如下:

其中:ε為非線性項(xiàng)中的小參數(shù),A為流體力參數(shù),它們可由實(shí)驗(yàn)數(shù)據(jù)獲得,一般情況下取ε=0.3,A=12,ωs為漩渦脫落頻率,St為Strouhal數(shù)。

1.3 立管和流體的耦合計(jì)算

將有內(nèi)流勻速流動(dòng)的頂張力立管的動(dòng)力學(xué)方程(1)和尾流振子模型方程(2)聯(lián)合,形成一個(gè)非線性振動(dòng)系統(tǒng)。采用Hermit插值函數(shù)將方程離散,得到立管振動(dòng)的矩陣方程形式:

本文采用Newmark-β方法,在時(shí)域里求解方程(3)、(4),可得立管的渦激振動(dòng)響應(yīng)。

為求解立管的橫向振動(dòng)動(dòng)力特性,將方程(3)忽略阻尼項(xiàng)及荷載項(xiàng)得

求解方程(5)可得立管渦激振動(dòng)的頻率和振型。

1.4 雨流計(jì)數(shù)法[10]

也叫塔頂計(jì)數(shù)法,是一種雙參數(shù)計(jì)數(shù)法,可以記錄載荷循環(huán)的全部信息,目前使用最為廣泛。其主要特點(diǎn)是根據(jù)材料的應(yīng)力—應(yīng)變過(guò)程進(jìn)行計(jì)數(shù),統(tǒng)計(jì)載荷波形中的循環(huán)和半循環(huán)。

根據(jù)計(jì)算得到的節(jié)點(diǎn)應(yīng)力—時(shí)程曲線,可獲得應(yīng)力—時(shí)間的峰谷值序列,本文采用4點(diǎn)法雨流計(jì)數(shù)原則進(jìn)行雨流計(jì)數(shù),計(jì)數(shù)條件如下:

(1)如果A>B,B≥D,C≤A;或者A

圖3 波形1Fig.3 Wave figure 1

(2)重復(fù)上述方法,取出一系列的全循環(huán),剩下的是發(fā)散—收斂序列,如圖5所示,已不再滿(mǎn)足上述計(jì)數(shù)條件,此時(shí)可用變程均值計(jì)數(shù)法,相鄰2個(gè)節(jié)點(diǎn)構(gòu)成一個(gè)半循環(huán)。范圍值Srange=|yi-yi+1|,幅值Sa=|yiyi+1|/2,均值Sm=|yi+yi+1|/2。再將具有相同均值和范圍的2個(gè)半波合成為一個(gè)全循環(huán)。

圖5 波形3Fig.5 Wave figure 3

(3)上述得到的所有全循環(huán)合在一起,完成對(duì)應(yīng)力—時(shí)程的雨流計(jì)數(shù)。

1.5 疲勞壽命分析

采用雨流計(jì)數(shù)法,得到應(yīng)力循環(huán)數(shù),然后運(yùn)用Palmgren-Miner線性累積疲勞損傷理論,計(jì)算立管的疲勞壽命。立管的疲勞累計(jì)損傷指數(shù)可以表示為:

其中:D是結(jié)構(gòu)的損傷指數(shù);ni是以Si為應(yīng)力幅值的實(shí)際循環(huán)次數(shù);N(si)是在應(yīng)力幅值為Si時(shí)的材料允許的極限循環(huán)次數(shù)。N(Si)可以通過(guò)S-N曲線取得,S-N曲線方程為:

其中:Δ σ為應(yīng)力幅值;loga為雙對(duì)數(shù)坐標(biāo)S-N曲線的截距;m為雙對(duì)數(shù)坐標(biāo)S-N曲線的斜率。

2 系統(tǒng)設(shè)計(jì)

系統(tǒng)以MATLAB 7.6.0為平臺(tái),使用圖形用戶(hù)界面開(kāi)發(fā)環(huán)境GUIDE,開(kāi)發(fā)出簡(jiǎn)潔清晰、使用方便的人機(jī)交互界面,以圖形、動(dòng)畫(huà)、數(shù)值等方式輸出分析結(jié)果。

圖形用戶(hù)界面(GUI)是使用圖形對(duì)象(如按鈕、文本框、滾動(dòng)條、菜單等)創(chuàng)建的用戶(hù)界面,用戶(hù)對(duì)每一個(gè)對(duì)象進(jìn)行編程,在GUI行為中達(dá)到相應(yīng)目的。實(shí)現(xiàn)一個(gè)GUI主要包括2項(xiàng)工作:GUI界面設(shè)計(jì)和GUI組件回調(diào)函數(shù)的編寫(xiě)[11]。

MA TLAB為用戶(hù)開(kāi)發(fā)圖形用戶(hù)界面提供了方便高效的集成環(huán)境:MATLAB圖形用戶(hù)界面開(kāi)發(fā)環(huán)境GUIDE(MATLAB’s Graphical User Iterface Development Environment)。圖形界面設(shè)計(jì)完成后,GUIDE可將其保存為2個(gè)文件:FIG文件和M文件[12]。

同一GUI不同callback間的數(shù)據(jù)傳遞有多種方法[13],如聲明全局變量、設(shè)置對(duì)象的userdata屬性、setappdata方法等。不同GUI之間的數(shù)據(jù)傳遞,可通過(guò)聲明全局變量實(shí)現(xiàn),本系統(tǒng)就是通過(guò)這種方式實(shí)現(xiàn)參數(shù)傳遞的,并可通過(guò)save、load等命令進(jìn)行保存和調(diào)用。

2.1 主要模塊

系統(tǒng)主要分為參數(shù)輸入模塊、分析計(jì)算模塊、結(jié)果輸出模塊,主界面見(jiàn)圖6。

圖6 主界面Fig.6 Main interface

圖7 立管參數(shù)Fig.7 Riser parameter

圖8 流體參數(shù)Fig.8 Fluid parameter

圖9 疲勞參數(shù)Fig.9 Fatigue parameter

2.1.1 參數(shù)輸入模塊 包括立管參數(shù)、流體參數(shù)和疲勞參數(shù)的輸入,界面見(jiàn)圖7～9所示。主要用到了Edit text,Static text,Pop-up Menu,ActiveX Control, Push Button等控件。系統(tǒng)通過(guò)聲明全局變量接收控件數(shù)據(jù),如通過(guò)全局變量riser_length獲取立管長(zhǎng)度:

riser_length=str2num(get(handles.riser_ length,’string’));

點(diǎn)擊確定時(shí)對(duì)數(shù)據(jù)合理性進(jìn)行識(shí)別,如判斷是否為數(shù)值、立管長(zhǎng)度是否等于水上長(zhǎng)度和水下長(zhǎng)度之和、外流速度截面深度是否大于水下長(zhǎng)度等,避免不必要的輸入錯(cuò)誤。

流體參數(shù)輸入界面中使用excel控件OWC11. Spreadsheet.11,用于接收外部流體速度數(shù)據(jù),并可查看流速截面,部分代碼如下:

spreadsheet=actxcontrol(’OWC11.Spreadsheet. 11’,[253 100 145 235]);%生成控件

current_profile=get(spreadsheet,’csvdata’);%獲取控件數(shù)據(jù),返回一個(gè)矩陣

參數(shù)輸入界面在打開(kāi)時(shí),系統(tǒng)對(duì)MATLAB工作空間進(jìn)行搜索,如變量已被導(dǎo)入工作空間,則將其值顯示在對(duì)應(yīng)的編輯框內(nèi),否則編輯框顯示為空,等待用戶(hù)輸入數(shù)據(jù)。如將工作空間中的全局變量riser_length值顯示在對(duì)應(yīng)編輯框內(nèi):

set(handles.riser_length,’string’,num2str(riser _length));

對(duì)于OWC11.Spreadsheet.11控件,MATLAB沒(méi)有有效的函數(shù)將數(shù)據(jù)填充在單元格內(nèi),本文自編了fill _spreadsheet函數(shù),實(shí)現(xiàn)了數(shù)據(jù)的直接顯示。

2.1.2 分析計(jì)算模塊 系統(tǒng)在計(jì)算過(guò)程中,對(duì)已知行列數(shù)的變量進(jìn)行預(yù)定義,提高了計(jì)算速度。通過(guò)編制myrainflow函數(shù),實(shí)現(xiàn)了雨流計(jì)數(shù)法在疲勞計(jì)算中的應(yīng)用。為使用戶(hù)更直觀地掌握計(jì)算進(jìn)度,對(duì)MATLAB提供的waitbar函數(shù)進(jìn)行了改進(jìn),編制mywaitbar函數(shù),在進(jìn)度條上顯示百分比。

2.1.3 結(jié)果輸出模塊

2.1.3.1 圖形結(jié)果 將立管任一點(diǎn)位移時(shí)程曲線、任一時(shí)刻位移曲線、最大位移包絡(luò)圖、任一點(diǎn)應(yīng)力時(shí)程曲線、最大應(yīng)力包絡(luò)圖、任一點(diǎn)功率譜圖、疲勞壽命集成在同一界面內(nèi),并設(shè)置互斥按鈕,可根據(jù)需要任意查看。

在任一時(shí)刻位移曲線選中狀態(tài)下,點(diǎn)擊動(dòng)畫(huà)演示可查看振動(dòng)動(dòng)畫(huà),輸入時(shí)刻為動(dòng)畫(huà)開(kāi)始時(shí)刻。各界面右下方,顯示立管的最大位移節(jié)點(diǎn)、最大應(yīng)力節(jié)點(diǎn)、最短疲勞壽命節(jié)點(diǎn)、最大位移、最大應(yīng)力和最短疲勞壽命,數(shù)據(jù)為只讀狀態(tài)。點(diǎn)擊圖形編輯進(jìn)入圖形編輯狀態(tài),可對(duì)坐標(biāo)和曲線任意編輯,也能保存成所需要的圖片格式。

2.1.3.2 數(shù)值結(jié)果 使用VideoSoft FlexArray Control控件顯示數(shù)值結(jié)果,包括立管各階固有頻率,各節(jié)點(diǎn)的最大位移、最大應(yīng)力和疲勞壽命。所有節(jié)點(diǎn)的最大位移、最大應(yīng)力和最短疲勞壽命將加粗顯示。數(shù)值結(jié)果只能查看,無(wú)法編輯。如:

set(h,’row’,i,’col’,1,’text’,num2str(natural_frequency(i)));%顯示節(jié)點(diǎn)的疲勞壽命

2.2 文件操作

系統(tǒng)將輸入?yún)?shù)和輸出結(jié)果保存到指定位置和名稱(chēng)的＊.mat數(shù)據(jù)庫(kù)文件中?？赏ㄟ^(guò)導(dǎo)入菜單,將其載入MATLAB工作空間,查看輸入?yún)?shù)和計(jì)算結(jié)果,也可修改參數(shù)并重新計(jì)算。

3 計(jì)算實(shí)例

本文選取的計(jì)算實(shí)例參數(shù)如下,立管兩端鉸接,總長(zhǎng)120 m,水下長(zhǎng)度100 m,水上長(zhǎng)度20 m,頂端張力200 kN,外徑0.25 m,內(nèi)徑0.211 6 m,立管密度7 700 kg/m3,立管材料彈性模量2.1e11N/m2,海水深度100 m,海水密度1 025 kg/m3,內(nèi)流密度800 kg/m3,附加質(zhì)量系數(shù)1,拖曳力系數(shù)0.7,其它參數(shù)如圖7～9所示。分別計(jì)算內(nèi)流速度為0,5,10,15,20 m/s時(shí)立管的動(dòng)力響應(yīng)及疲勞壽命,結(jié)果見(jiàn)表1。圖10～17列出了內(nèi)流速度為0 m/s時(shí)系統(tǒng)的計(jì)算結(jié)果截圖。

表1 不同內(nèi)流速度下計(jì)算結(jié)果Table 1 Results of different internal flow velocity

圖10 任一點(diǎn)位移時(shí)程曲線Fig.10 Displacement-time curve

圖11 任一時(shí)刻位移曲線Fig.11 Displacement curve

圖12 最大位移包絡(luò)圖Fig.12 Envelope of displacement

圖13 任一點(diǎn)應(yīng)力時(shí)程曲線Fig.13 Stress-time curve

圖14 最大應(yīng)力包絡(luò)圖Fig.14 Envelope of stress

圖15 任一點(diǎn)功率譜圖Fig.15 Power spectrum

圖16 疲勞壽命Fig.16 Fatigue life

圖17 數(shù)值結(jié)果Fig.17 Numerical result

計(jì)算結(jié)果表明,隨著內(nèi)流速度的增加,立管固有頻率降低。隨著立管固有頻率的變化,當(dāng)其接近漩渦脫落頻率時(shí),立管將進(jìn)入鎖振區(qū)域,增大渦激振動(dòng)振幅,縮短疲勞壽命。同樣,立管固有頻率的變化也可使立管脫離鎖振區(qū)域,減小渦激振動(dòng)振幅,延長(zhǎng)疲勞壽命。因此,內(nèi)流對(duì)立管的渦激振動(dòng)響應(yīng)和疲勞壽命均有重要影響,在工程設(shè)計(jì)中應(yīng)予以考慮。

4 結(jié)語(yǔ)

本文考慮內(nèi)流作用,基于MATLAB開(kāi)發(fā)出《海洋立管渦激振動(dòng)數(shù)值模擬系統(tǒng)NSVIV 1.0》,分析立管在海流作用下的渦激振動(dòng)響應(yīng),得到立管的橫向振動(dòng)特性及疲勞壽命。系統(tǒng)界面簡(jiǎn)潔清晰、使用方便,適用于頂張力立管的分析,隨著研究的不斷深入,將進(jìn)一步升級(jí)為功能齊全的海洋立管渦激振動(dòng)分析軟件,最終有望用于我國(guó)海洋立管的工程設(shè)計(jì)。

[1] 秦延龍,王世澎.海洋立管渦激振動(dòng)計(jì)算方法進(jìn)展[J].中國(guó)海洋平臺(tái).2008,23(4):14-17.

[2] J.Kim Vandiver,Li Li.SHEAR7 V4.4 PROGRAM THEORETICAL MANUAL[R].Cambridge,MIT:2005.

[3] Michael Triantafyllou.VIVA MARINE RISER SOFTWARE USERS MANUAL WINDOWS VERSION 2.0[R].Cambridge: MIT,2000.

[4] Larsen C M,Vikestad K,Yttevik R,et al.VIVANA,THEORY MANUAL[R].Norway:Trondheim,2001.

[5] 郭海燕,王樹(shù)青,劉德輔.海洋環(huán)境荷載下輸液立管的靜、動(dòng)力特性研究[J].青島海洋大學(xué)學(xué)報(bào),2001,31(4):605-611.

[6] GUO Hai-yan,LOU Min.Effect of internal flow on vortex-induced vibration of risers[J].Journal of Fluids and Structures, 2008,24(4):496-504.

[7] 郭海燕,董文乙,婁敏.海中輸流立管渦激振動(dòng)試驗(yàn)研究及疲勞壽命分析[J].中國(guó)海洋大學(xué)學(xué)報(bào):自然科學(xué)版,2008,38(3): 503-507.

[8] Guo Hai-yan,Li Xiao-min,Liu Xiao-chun.Numerical prediction of vortex induced vibrations on top tensioned riser in consideration of internal flow[J].China Ocean Engineering,2008,22(4):675-682.

[9] Facchinetti M L,Langre E D,Biolley F.Coupling of structure and wake oscillators in vortex-induced vibrations[J].Journal of Fluids and Structures,2004,19(2):123-140.

[10] Anzai H,Endo T.On-site fatigue damage under complex loading [J].International Journal of Fatigue,1979,1(1),49-57.

[11] 施曉紅,周佳.精通GUI圖形界面編程[M].北京:北京大學(xué)出版社,2003.

[12] 陳杰.MATLAB寶典[M].北京:電子工業(yè)出版社,2008.

[13] 黃菲,劉振興,羅銘.基于圖形用戶(hù)界面(GUI)的異步電動(dòng)機(jī)仿真系統(tǒng)[J].虛擬儀器與應(yīng)用,2008,31(6):178-180.

Abstract: Numerical Simulation System of Vortex-Induced Vibration of Marine Riser(NSVIV 1.0) based on MATLAB is developed.Wake oscillator model is used to simulate the effect of external fluid on the riser and the effect of the internal flowing fluid on the riser is considered.VIV dynamic responses and fatigue life of riser is predicted and analyzed.The interface of system is concise,clear,easy to use and suitable for TTR,and lays the foundation for developing more full-featured design and analysis software for marine riser.

Key words: MATLAB;GUI;marine riser;Vortex-Induced Vibration(VIV);numerical simulation

責(zé)任編輯 陳呈超

Research on Numerical Simulation System of Vortex-Induced Vibration of Marine Riser Based on MATLAB

LI Hong-Chun,GUO Hai-Yan,LI Xiao-Min
(College of Engineering,Ocean University of China,Qingdao 266100,China)

TV312

A

1672-5174(2010)09Ⅱ-207-06

國(guó)家高技術(shù)研究發(fā)展計(jì)劃項(xiàng)目(2006AA09Z356,2007AA09Z313)資助

2009-07-14;

2009-10-30

李洪春(1983-),男,碩士生。E-mail:lhc627@163.com