蔺云飞,肖文生,崔俊国,马一心,施海滨,苑思敏
(1. 中国石油大学(华东) 机电工程学院,山东 青岛 266555;2. 上海振华重工(集团)股份有限公司,上海 200125)
风能作为一种可再生能源,近年来得到快速的发展,其中海上风电的发展也取得了巨大的进步[1]。目前世界各国对沿海风电,尤其是深远海风能资源开发都尤为重视[2],虽然深海区域的风能资源充沛,但近年来海上风电的发展多在近海和浅海(水深不超过50 m),其主要原因是深海情况下风机组施工难度大,设备运行环境恶劣,现有的一些风电安装平台不能适应深水作业,为此创新设计一种基于正压型下浮体的3500 型深水高效海上风电安装平台,如图1所示。
图 1 3500 型风电安装平台结构原理图Fig. 1 Structure schematic diagram of type 3500 wind power installation platform
3500 型海上风电安装平台通过底部中空的正压型下浮体冲排水获得浮力顶升平台实现作业,作业期间,承受浮体内超压和外界水压。下浮体是一个体积超20000 m3的密闭区域,且表面连接有输气管线等部件,而在实际工程中容器气密性越高,则建造难度和成本越高,为此下浮体允许制造过程有微量的泄漏,但泄漏过程要满足下浮体建立保压控制系统的要求[3]。要想实现下浮体加压及保压控制系统的设计,对下浮体泄漏过程进行研究,即控制下浮体的气密性非常有必要。
目前,对于气密性的研究多为气密性检测方法,如气泡法和涂抹法等[4],且对气密性泄漏过程压降规律的研究更多是对泄漏量和泄漏压力的理论计算,而对泄漏过程影响因素的探究相对较少。黄震等[5]通过理论计算方法研究了航天领域中密封舱再入过程中的泄漏压降规律。王兆芹[6]分析了高压输气管道泄漏孔径与管径之比对泄漏速率的影响。Schatzel 等[7]通过实验研究了煤矿井封件的泄漏速率。本文建立下浮体容器泄漏过程的数学模型,结合Fluent 对下浮体理论计算结果进行仿真分析验证,并探究下浮体泄漏口位置和数量以及下浮体内部结构对压降过程的影响。
下浮体计算时不考虑箱体变形所带来的容积变化,即下浮体静态泄漏模型可简化为图2 所示的结构,模型容积为V,内压力为Pi,气体密度为ρi,外界环境压力为Pe,密度为ρe,假设所有泄漏口都由一个面积为A的等效泄漏孔来代替,开口处压力和密度分别为P2和ρ2,流出速度为u2。
图 2 泄漏模型图Fig. 2 Leakage model diagram
下浮体在发生泄漏时,可以将泄漏过程看成许多微小时间段dt组合而成,而每个时间段dt都处于稳态过程,并且由于泄漏口尺寸相对下浮体的尺寸较小,整个过程可以看成可逆绝热过程,为此通过以下流体运动基本方程推导下浮体泄漏过程的压降参数。
式中:z为流道壁对流体得切向摩擦应力;L为流道周长;δq为流体与外界交换的热量;G为重力加速度;δwnen为外界做的功。
基于等熵可压缩流动理论的压降参数模型为[8]:
基于不可压缩流动理论的压降参数模型为:
式中,Cp为形状系数。
本文研究的3500 型风电安装平台主要在我国南海沿岸区域工作,以往风电安装设备多是针对浅海区域,业内常把水深60 m 以内作为浅海区域,而此平台可适应深水60 m 及以上区域。下浮体工作时要求内部保持正压,即保证内压Pi大于外压Pe,压差不超过一个大气压,研究时下浮体泄漏口取光滑圆孔,此时可以得到下浮体初始状态时的相关参数,如表1 所示。
表 1 初始泄漏参数Tab. 1 Initial leak parameter
结合表1 数据,可以计算泄漏开始时的流速和马赫数值:
式中:Ma为马赫数;c为声速,在水中取1500 m/s。
一般来说没有不可压缩的流体,经相关计算表明,当0
通过计算可以知道下浮体泄漏过程中马赫数小于0.3,故可忽略其压缩性,计算采用不可压缩流动理论,对式(6)解微分方程并代入t=0 时压差∆Pman得:∆p
式中,为t时刻下浮体内外压差。
为研究下浮体装置的不同泄漏口面积和水深的压降情况,计算下浮体达到平衡状态所需时间,平衡状态即压差时 ∆P=0,即可以得到平衡时间t与水深h和泄漏面积A的函数关系式:
为更好地分析下浮体泄漏过程压降规律,根据下浮体的水深作业范围,探究下浮体在60 m,70 m,80 m,90 m 和100 m 等5 类水深下平衡时间与泄漏面积的关系,如图3 所示。
图 3 平衡时间与泄漏面积关系Fig. 3 Relationship between equilibrium time and leakage area
可以看出,在相同水深情况下随着泄漏面积的增加,下浮体内压到达平衡状体的时间逐渐增大,而时间增大的速度却随着泄漏面积的增多而逐渐减少。此外,下浮体泄漏面积小于0.1 m2时,所需平衡时间趋于无穷大,而泄漏面积大于1 m2时,所需平衡时间基本不再受泄漏面积影响,趋于一个稳定值。为此对下浮体气密性泄漏面积的研究主要集中在0.1~1 m2,即研究0.1 m2,0.25 m2,0.5 m2,0.75 m2和1 m2等5 五种泄漏面积下的压降规律,表2 为不同状态下下浮体到达平衡状态所需时间。
表 2 不同状态下达到平衡所需时间Tab. 2 Time required to reach equilibrium in different states
采用仿真模拟方法进行数学模型的验证,试验条件基于工况和正压型下浮体的特性而定,而下浮体泄漏面积需要制作等尺寸模型进行试验测得,本次验证以0.5 m2的圆形泄漏口和作业水深60 m 为试验条件,下浮体仿真模型长宽高分别为90 m,48 m 和5 m,并开有4 个半径5 m 的圆柱桩腿孔,如图4 所示。
图 4 仿真模型结构图Fig. 4 Structure diagram of simulation model
计算采用Ansys Fluent 模块进行数值模拟计算,仿真过程中,为保证下浮体泄漏口处的压力在泄漏过程中保持不变,设置泄漏口出口压力为0 Pa,操作压力取工作水深的外压,初始化后通过patch 给内部计算区域施加初始压力101325 Pa。为直观显示下浮体内部计算区域的压力分布情况,在下浮体内部沿3 个方向建立观测截面,图5 为各观察截面的压力云图。
图 5 各截面压力云图Fig. 5 Pressure cloud of each section
可以看出,靠近泄漏口区域处压力变化较大,与周边区域压力差超过1000 Pa,此时应避免在泄漏口附近建立监测点,转而在其周围区域建立监测点。图5(b)可以清楚表现出其他区域处压力具体范围,可以看出除泄漏口区域,各单元之间压力变化小,压差不超过1Pa,为此可以在下浮体泄漏口区域外的任意区域取一点作为压力监测点,通过监测点的压力变化来反映下浮体泄漏过程中的整体压力变化。
为了减小网格尺寸引起的计算误差并选取合理的网格尺寸,对计算模型进行网格独立性检验[10],整个计算区域采用四面体进行网格划分并添加膨胀层,在泄漏口部分区域进行网格加密,如图6 所示。
图 6 泄漏口处网格划分Fig. 6 Mesh division at leakage outlet
表3 给出了用于计算模型内部压力变化情况的4 种不同网格参数,编号从Mesh1~Mesh4,4 种网格均划分了膨胀层,最大网格尺寸取0.4 m,0.5 m 和0.6 m,膨胀层分为4 层和5 层进行网格的验证,网格数量从107 万增加到292 万,网格质量以Skewness 值衡量,高Skewness 的单元极易造成计算的发散,而其值越接近0,则说明网格质量越好。
表 3 网格尺寸Tab. 3 Size of mesh opening
图7 表示4 种不同网格尺寸下,下浮体经Fluent 仿真模拟后内部压力变化曲线。可以看出,下浮体达到平衡状态的时间从Mesh1~Mesh4呈增大趋势,Mesh1 到达平衡时间为430.5 s,Mesh2 为438.0 s,Mesh3 为444.0 s,Mesh4 为446.5 s,而理论公式计算平衡时间为434.1 s。
图 7 不同网格尺寸下的压降曲线Fig. 7 Pressure drop curves under different mesh sizes
为更直观地反映各曲线与理论公式的相似性,通过Spss 软件进行数据分析计算数据间相关性,相关性计算值为Pearson,Spearman 和Kendall 系数,如表4 所示。3 个系数反映的都是2 个变量之间变化趋势的方向以及程度,其绝对值越接近1,相关性就越强,越接近0,则相关性越弱[11]。
表 4 各曲线间相关性系数Tab. 4 Correlation Coefficient between curves
在仿真分析过程中Mesh1 平衡时间与Mesh2 平衡时间的增量为1.86%,Mesh2 与Mesh3 增量为1.37%,而Mesh3 加密至Mesh4 时增量仅变化了0.4%。可以看出,Mesh3 和Mesh4 的相关系数基本一致。仿真计算结果和理论公式计算结果较为吻合,验证了理论公式与仿真分析的正确性,同时也可以确定Mesh3 为网格尺寸的基准进行后续分析。
对泄漏口数量进行计算时,模型泄漏口包含1~4个的情况,面积关系为S1=2S2=3S3=4S4=0.5 m2;泄漏口位置分析时分别在上表面(above)、前表面(front)及左表面(left)进行开口,同时为了便于观察和记录,对泄漏口数量和位置进行编号,如上表面开1 个泄漏口表示为above-1。图8 为下浮体仿真模型各表面不同泄漏口数量下的压降曲线。
图 8 各表面不同泄漏口数压降曲线Fig. 8 Pressure drop curves of different leakage ports on the same plane
可以判断,随着泄漏口数的增多,下浮体内部到达平衡的时间逐渐减少,为探究这一现象取y=0 这一截面分析其流速图,如图9 所示。可以看出,不同数量泄漏口时下浮体内分为不同的气体区域,泄漏口数越多气体区域也越多,相校于泄漏口数少的模型,气体分子流经的区域也就越小,所以流出的时间也会稍微减少。
图 9 不同泄漏口数流速图Fig. 9 Flow rate diagram of different leakage ports
为更好分析不同压降曲线的差异程度,引入均方根误差(RMSE)指标,计算时以上一曲线为基准,计算公式为:
通过计算可以知道各表面不同泄漏口数下平衡时间的时差及数据间的RMSE 值,如表5 所示。表中平衡时间为整个泄漏过程结束时刻,如果把平衡时间取为433.5s,那么1 个和2 个泄漏口之间时差占总时间的1.2%~1.3%,此时RMSE值为450 左右;2 和3 个泄漏口之间时差占0.90%~0.97%,此时RMSE值为230左右;3 和4 个泄漏口之间时差仅占0.39%~0.43%,此时RMSE值不超过100。
表 5 各表面不同泄漏口数下时差和RMSETab. 5 The time difference and rmse of different leakage ports on each surface
即在同一表面内随着泄漏口数的增多,平衡时差和RMSE值都逐渐减小且减小幅度变缓,而3 个和4 个泄漏口时的压降情况基本一致。可以认为当泄漏口数大于一定量时,压力变化曲线与泄漏口数量无关,同时可以把RMSE值作为压降过程差异大小的标准,当2 个过程RMSE小于100 时,压降情况可以视为一致。
对不同表面与不同泄漏口数情况下的压降规律进行差异分析,从图10 可以看出,泄漏口数相同情况下不同表面间的曲线存在些许差异。
图 10 各泄漏口数不同表面压降曲线Fig. 10 Different surface pressure drop curves for each leakage port number
为进一步验证差异,计算各表面在泄漏口数相同时压降数据的RMSE值,如表6 所示。其中不同表面间相同泄漏口数下的时差最大为1.8 s,占平衡时间的0.41%。同时各泄漏过程间RMSE均小于100,此时就可以把各压降情况视为一致,即泄漏口位置对泄漏过程的影响几乎可以忽略不计。
表 6 各状态间的平衡时差和RSME 值Tab. 6 Equilibrium time difference and rsme value between states
下浮体上下箱体中间有横向或纵向加强筋板支撑结构,图11 为下浮体内部筋板结构。为探讨下浮体内部结构对泄漏过程的影响,对比下浮体内部无筋板,横向筋板和纵向筋板3 种情况,其中在不同下浮体内部结构中横向筋板和纵向筋板数量一致,筋板体积占下浮体总体积的0.3%,可以忽略不计。
下浮体内部无筋板,有横向筋板和纵向筋板的压降曲线如图12 所示。可以看出在有筋板的情况下横向筋板压降和纵向筋板压降情况相差较小,无筋板和有筋板情况相差也不大。
图 11 内部筋板结构Fig. 11 Internal stiffened plate construction
图 12 不同内部结构下的影响曲线Fig. 12 Influence curves of different internal structures
表7 为3 种内部结构压降平衡时间的时差和压降过程的RMSE值。可以看出无筋板的情况下与横向筋板和纵向筋板的时差占比0.39%和0.48%,横向筋板与纵向筋板的时差占比仅为0.14%,并且RMSE值均小于100。即在筋板体积占比下浮体总体积较小的情况下,无筋板泄漏,横向筋板泄漏和纵向筋板泄漏的泄漏状态基本保持一致。
表 7 不同内部结构的平衡时差和RMSE 值Tab. 7 Balance time difference and rmse values for different internal structures
基于Fluent 仿真模拟,将有限元仿真结果与数学模型理论计算结果对比可知结果较为吻合。对下浮体泄漏过程中不同影响因素进行研究,结果表明:随着泄漏口数量的增多,下浮体平衡时间逐渐减少,并且3 个和4 个泄漏口的压降过程基本保持一致,即泄漏口数量超过一定量时压降过程不会发生变化;泄漏口所在的表面位置不影响泄漏过程,且在下浮体内部结构体积占比较小情况下,内部筋板结构对压降过程影响可以忽略不计。
本文研究为下浮体密封和加压及保压控制系统设计提供了理论基础,同时也为相关领域探究容器泄漏过程提供了参考依据。