基于MA TLAB的海洋立管涡激振动数值模拟系统研究

2010-09-13李洪春郭海燕李效民

中国海洋大学学报（自然科学版） 2010年1期

关键词：涡激立管控件

李洪春,郭海燕,李效民

(中国海洋大学工程学院,山东青岛266100)

基于MA TLAB的海洋立管涡激振动数值模拟系统研究

李洪春,郭海燕,李效民

(中国海洋大学工程学院,山东青岛266100)

基于MATLAB开发海洋立管涡激振动数值模拟系统NSVIV 1.0,系统采用尾流振子模型模拟外流对立管结构的作用,考虑内流对立管结构的影响,对立管的涡激振动动力响应和疲劳寿命进行预测分析。系统界面简洁清晰,使用方便,适用于顶张力立管,为进一步集成功能齐全的海洋立管设计分析软件打下基础。

MATLAB;图形用户界面;海洋立管;涡激振动;数值模拟

立管是海面作业平台与海底井口间的主要连接件,是海洋油气开发的重要设施。立管运行期间受力复杂,当波浪或海流流经立管时,在一定的流速条件下,可在立管两侧交替形成旋涡,当旋涡脱落频率与立管固有频率接近时,将使立管强烈振动,发生“锁定”(lock-in)现象,加剧立管的疲劳破坏。因此,涡激振动引起的疲劳损伤,是立管设计中必须考虑的重要因素。

立管涡激振动的数值模拟研究,主要有2种方法[1]:一种是通过求解Navier-Strokes方程的计算流体动力学法(CFD)。另一种是通过实验数据确定流体力参数的半经验法,基于此种方法的几种计算模型在海洋工程中得到广泛的应用,比较著名的有MIT的Vandiver教授开发的SHEAR7[2]、MIT的Triantafyllou教授开发的VIVA[3]、挪威的Larsen等开发的VIVANA[4]等。

内流对立管振动的影响已引起愈来愈多学者的重视,郭海燕等[5-7]考虑管内流体和管外海洋环境荷载的共同作用,分析立管的静、动力特性,并开展了相关理论和实验研究。

目前国外的涡激振动分析软件未考虑立管内部流体的作用,本文开发的海洋立管涡激振动数值模拟系统—NSVIV 1.0,考虑内流作用,采用半经验的尾流振子模型模拟流体力,用时程分析法,求解立管在任意剖面形式海流作用下的涡激振动响应和疲劳寿命,计算结果将更接近工程实际。系统基于常用数学软件MATLAB开发,界面简洁清晰,使用方便。NSVIV 1.0适用于顶张力立管的分析,随着研究的不断深入,将升级为包括钢悬链线立管等在内、功能齐全的海洋立管涡激振动分析软件,最终有望用于我国海洋立管的工程设计。

1 计算原理

1.1 立管动力学模型

立管受内部流体和外部海洋环境荷载的共同作用,见图1。假定管外海流沿X轴正向流动,管内流体自上而下流动。立管将在YOZ平面内产生横向涡激振动,在XOZ平面内产生顺流向振动,本系统主要研究立管横向的动力响应,取立管单元dz进行分析,见图2。

图1 海洋立管振动示意图Fig.1 Vibration of marine riser

图2 立管单元dz变形图Fig.2 Deformation of riser element dz

根据功-能平衡原理,可得立管横向涡激振动方程[8]:

其中:m=mr+mi+ma,mr为单位长度立管质量,mi为单位长度管内流体质量,ma为附加流体质量;y为立管横向位移;c为结构阻尼系数,c′为流体阻尼系数;V为内部流体流速;E为立管弹性模量,I为立管截面惯性矩;Te表示立管单元的有效轴向力;ρsea为海水密度,U为外流流速;D为立管直径;CL0为升力系数幅值,一般取CL0=0.3;q为升力振子系数,q=2CL/CL0。

1.2 流体力模型

本文采用Facchinetti et al.[9]推荐的尾流振子模型,该模型只考虑流体与结构之间的线性耦合项,其表达式如下:

其中:ε为非线性项中的小参数,A为流体力参数,它们可由实验数据获得,一般情况下取ε=0.3,A=12,ωs为漩涡脱落频率,St为Strouhal数。

1.3 立管和流体的耦合计算

将有内流匀速流动的顶张力立管的动力学方程(1)和尾流振子模型方程(2)联合,形成一个非线性振动系统。采用Hermit插值函数将方程离散,得到立管振动的矩阵方程形式:

本文采用Newmark-β方法,在时域里求解方程(3)、(4),可得立管的涡激振动响应。

为求解立管的横向振动动力特性,将方程(3)忽略阻尼项及荷载项得

求解方程(5)可得立管涡激振动的频率和振型。

1.4 雨流计数法[10]

也叫塔顶计数法,是一种双参数计数法,可以记录载荷循环的全部信息,目前使用最为广泛。其主要特点是根据材料的应力—应变过程进行计数,统计载荷波形中的循环和半循环。

根据计算得到的节点应力—时程曲线,可获得应力—时间的峰谷值序列,本文采用4点法雨流计数原则进行雨流计数,计数条件如下:

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

图3 波形1Fig.3 Wave figure 1

图5 波形3Fig.5 Wave figure 3

(3)上述得到的所有全循环合在一起,完成对应力—时程的雨流计数。

1.5 疲劳寿命分析

采用雨流计数法,得到应力循环数,然后运用Palmgren-Miner线性累积疲劳损伤理论,计算立管的疲劳寿命。立管的疲劳累计损伤指数可以表示为:

其中:D是结构的损伤指数;ni是以Si为应力幅值的实际循环次数;N(si)是在应力幅值为Si时的材料允许的极限循环次数。N(Si)可以通过S-N曲线取得,S-N曲线方程为:

其中:Δ σ为应力幅值;loga为双对数坐标S-N曲线的截距;m为双对数坐标S-N曲线的斜率。

2 系统设计

系统以MATLAB 7.6.0为平台,使用图形用户界面开发环境GUIDE,开发出简洁清晰、使用方便的人机交互界面,以图形、动画、数值等方式输出分析结果。

图形用户界面(GUI)是使用图形对象(如按钮、文本框、滚动条、菜单等)创建的用户界面,用户对每一个对象进行编程,在GUI行为中达到相应目的。实现一个GUI主要包括2项工作:GUI界面设计和GUI组件回调函数的编写[11]。

MA TLAB为用户开发图形用户界面提供了方便高效的集成环境:MATLAB图形用户界面开发环境GUIDE(MATLAB’s Graphical User Iterface Development Environment)。图形界面设计完成后,GUIDE可将其保存为2个文件:FIG文件和M文件[12]。

同一GUI不同callback间的数据传递有多种方法[13],如声明全局变量、设置对象的userdata属性、setappdata方法等。不同GUI之间的数据传递,可通过声明全局变量实现,本系统就是通过这种方式实现参数传递的,并可通过save、load等命令进行保存和调用。

2.1 主要模块

系统主要分为参数输入模块、分析计算模块、结果输出模块,主界面见图6。

图6 主界面Fig.6 Main interface

图7 立管参数Fig.7 Riser parameter

图8 流体参数Fig.8 Fluid parameter

图9 疲劳参数Fig.9 Fatigue parameter

2.1.1 参数输入模块包括立管参数、流体参数和疲劳参数的输入,界面见图7～9所示。主要用到了Edit text,Static text,Pop-up Menu,ActiveX Control, Push Button等控件。系统通过声明全局变量接收控件数据,如通过全局变量riser_length获取立管长度:

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

点击确定时对数据合理性进行识别,如判断是否为数值、立管长度是否等于水上长度和水下长度之和、外流速度截面深度是否大于水下长度等,避免不必要的输入错误。

流体参数输入界面中使用excel控件OWC11. Spreadsheet.11,用于接收外部流体速度数据,并可查看流速截面,部分代码如下:

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

current_profile=get(spreadsheet,’csvdata’);%获取控件数据,返回一个矩阵

参数输入界面在打开时,系统对MATLAB工作空间进行搜索,如变量已被导入工作空间,则将其值显示在对应的编辑框内,否则编辑框显示为空,等待用户输入数据。如将工作空间中的全局变量riser_length值显示在对应编辑框内:

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

对于OWC11.Spreadsheet.11控件,MATLAB没有有效的函数将数据填充在单元格内,本文自编了fill _spreadsheet函数,实现了数据的直接显示。

2.1.2 分析计算模块系统在计算过程中,对已知行列数的变量进行预定义,提高了计算速度。通过编制myrainflow函数,实现了雨流计数法在疲劳计算中的应用。为使用户更直观地掌握计算进度,对MATLAB提供的waitbar函数进行了改进,编制mywaitbar函数,在进度条上显示百分比。

2.1.3 结果输出模块

2.1.3.1 图形结果将立管任一点位移时程曲线、任一时刻位移曲线、最大位移包络图、任一点应力时程曲线、最大应力包络图、任一点功率谱图、疲劳寿命集成在同一界面内,并设置互斥按钮,可根据需要任意查看。

在任一时刻位移曲线选中状态下,点击动画演示可查看振动动画,输入时刻为动画开始时刻。各界面右下方,显示立管的最大位移节点、最大应力节点、最短疲劳寿命节点、最大位移、最大应力和最短疲劳寿命,数据为只读状态。点击图形编辑进入图形编辑状态,可对坐标和曲线任意编辑,也能保存成所需要的图片格式。

2.1.3.2 数值结果使用VideoSoft FlexArray Control控件显示数值结果,包括立管各阶固有频率,各节点的最大位移、最大应力和疲劳寿命。所有节点的最大位移、最大应力和最短疲劳寿命将加粗显示。数值结果只能查看,无法编辑。如:

set(h,’row’,i,’col’,1,’text’,num2str(natural_frequency(i)));%显示节点的疲劳寿命

2.2 文件操作

系统将输入参数和输出结果保存到指定位置和名称的＊.mat数据库文件中。可通过导入菜单,将其载入MATLAB工作空间,查看输入参数和计算结果,也可修改参数并重新计算。

3 计算实例

本文选取的计算实例参数如下,立管两端铰接,总长120 m,水下长度100 m,水上长度20 m,顶端张力200 kN,外径0.25 m,内径0.211 6 m,立管密度7 700 kg/m3,立管材料弹性模量2.1e11N/m2,海水深度100 m,海水密度1 025 kg/m3,内流密度800 kg/m3,附加质量系数1,拖曳力系数0.7,其它参数如图7～9所示。分别计算内流速度为0,5,10,15,20 m/s时立管的动力响应及疲劳寿命,结果见表1。图10～17列出了内流速度为0 m/s时系统的计算结果截图。

表1 不同内流速度下计算结果Table 1 Results of different internal flow velocity

图10 任一点位移时程曲线Fig.10 Displacement-time curve

图11 任一时刻位移曲线Fig.11 Displacement curve

图12 最大位移包络图Fig.12 Envelope of displacement

图13 任一点应力时程曲线Fig.13 Stress-time curve

图14 最大应力包络图Fig.14 Envelope of stress

图15 任一点功率谱图Fig.15 Power spectrum

图16 疲劳寿命Fig.16 Fatigue life

图17 数值结果Fig.17 Numerical result

计算结果表明,随着内流速度的增加,立管固有频率降低。随着立管固有频率的变化,当其接近漩涡脱落频率时,立管将进入锁振区域,增大涡激振动振幅,缩短疲劳寿命。同样,立管固有频率的变化也可使立管脱离锁振区域,减小涡激振动振幅,延长疲劳寿命。因此,内流对立管的涡激振动响应和疲劳寿命均有重要影响,在工程设计中应予以考虑。

4 结语

本文考虑内流作用,基于MATLAB开发出《海洋立管涡激振动数值模拟系统NSVIV 1.0》,分析立管在海流作用下的涡激振动响应,得到立管的横向振动特性及疲劳寿命。系统界面简洁清晰、使用方便,适用于顶张力立管的分析,随着研究的不断深入,将进一步升级为功能齐全的海洋立管涡激振动分析软件,最终有望用于我国海洋立管的工程设计。

[1] 秦延龙,王世澎.海洋立管涡激振动计算方法进展[J].中国海洋平台.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] 郭海燕,王树青,刘德辅.海洋环境荷载下输液立管的静、动力特性研究[J].青岛海洋大学学报,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] 郭海燕,董文乙,娄敏.海中输流立管涡激振动试验研究及疲劳寿命分析[J].中国海洋大学学报:自然科学版,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].北京:北京大学出版社,2003.

[12] 陈杰.MATLAB宝典[M].北京:电子工业出版社,2008.

[13] 黄菲,刘振兴,罗铭.基于图形用户界面(GUI)的异步电动机仿真系统[J].虚拟仪器与应用,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

责任编辑陈呈超

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

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

国家高技术研究发展计划项目(2006AA09Z356,2007AA09Z313)资助

2009-07-14;

2009-10-30

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