基于精细化建模的换热管涡激振动数值研究

孙军帅, 李旭鹏, 孙汝雷, 温济铭, 李小畅, 田瑞峰

(1.哈尔滨工程大学 核科学与技术学院, 黑龙江 哈尔滨 150001; 2.中国船舶及海洋工程设计研究院, 上海 200011; 3.哈尔滨工程大学 航天与建筑工程学院, 黑龙江 哈尔滨 150001)

核电站蒸汽发生器作为核电站关键设备,其安全性关乎核电站的安全运行。在长期运行中,蒸汽发生器换热管极容易受流体作用而发生振动,当振动剧烈时,换热管会与相邻管束、支承结构等发生碰撞或摩擦,导致管束失效,造成巨大经济损失[1],有必要对换热管的流致振动现象进行深入研究。

大涡模拟方法(large eddy simulation, LES)在湍流计算方面具有优势,随着计算机技术的迅猛发展,具有较好的使用前景[2]。Lourenco等[3]较早地采用了粒子图像测速法(particle image velocimetry, PIV)技术展开了Re=3 900的圆柱绕流研究。文献[4-5]采用PIV技术对尾流区流场以及剪切层不稳定性等进行了研究。由于流场的湍流特性受实验条件影响较大,且通过实验研究能够捕获的流场信息极为有限,因此大量学者通过数值手段展开相关研究。Beaudan等[6]较早地采用LES方法对Re=3 900的圆柱绕流问题进行了数值模拟。Kravchenko等[7]基于B样条方法的高精度LES方法,研究了Re=3 900圆柱绕流问题,通过与实验结果[3]对比发现,实验可能存在环境干扰从而导致了剪切层较早转捩。Parnaudeau等[5]除了进行了实验研究外,同样采用LES方法展开了数值研究。战庆亮等[8]、端木玉等[9]采用LES方法研究了Re=3 900时圆柱绕流流场特性,通过对比研究验证了计算程序在湍流场计算的准确性。一些关键流场参数方面仍存在一定的争议,这可能是由于实验技术差异以及数值计算模型不同导致的。

单根换热管流致振动从流致振动机理来说,本质上属于单圆柱涡激振动的范畴。文献[10-13]对圆柱涡激振动进行了系统性的实验研究。文献[14-15]在研究低质量阻尼参数的双自由度圆柱涡激振动时发现:当质量比m*=2.6时,横向振动响应幅值会产生“超级上端分支”,并出现了“2 T”形态的尾迹模式。Kang等[16]、Li等[17]基于RANS法建立了双自由度圆柱涡激振动问题的二维计算模型,但由于二维模型忽视了流场在高雷诺数时具有高度三维化涡体结构,导致得到的最大横向振幅和“超级上端分支”的范围与实验结果存在一定的偏差。Gsell等[18]和Li等[17]分别采用直接数值模拟(direct numerical simulation,DNS)方法和LES方法建立了双自由度圆柱涡激振动问题的三维计算模型以模拟“超级上端分支”,但结果还是均低于实验值;同时,由于计算工况较少,无法对该问题进行较好描述。

综上所述,在单圆柱涡激振动问题上,主要存在以下问题:1)在流场计算方面,CFD模型存在差异,导致关键流场参数存在着一定争议;2)在振动响应计算方面,不论建立的计算模型是二维还是三维的,计算结果均与Jauvtis &Williamson[15]的实验值存在较大差异,体现在超级上端分支范围、最大横向振幅、运动轨迹等方面。

本文基于LES方法建立了精细化流场计算模型,以雷诺数为3 900的固定圆柱绕流为例,进行了网格敏感性分析,而后通过不同求解方法的对比研究选取HHT-α方法求解结构动力学方程,最终结合局部动网格变形技术,建立了换热管涡激振动流固耦合计算模型,并对Jauvtis &Williamson[15](简称为J-W)的经典涡激振动实验的振动响应进行了预测。

1 流固耦合计算模型

1.1 大涡模拟控制方程

LES方法最早是由气象学家Smagorinsky[19]提出的,其主要思路是对N-S方程进行滤波处理,将湍流中涡分为大尺度涡和小尺度涡2种,对大涡直接计算,对小涡进行模拟。对于不可压缩流动,滤波后的N-S方程为:

(1)

(2)

亚格子应力模型有很多形式[19-21],本文采用了WALE模型[21],该模型通过考虑湍流壁面及动量传递效应等,保证了对近壁面流域模拟的准确性。同时,采用SIMPLEC方法求解压力速度耦合方程,采用二阶隐式方法对时间项进行离散求解。

1.2 结构动力学模型

将单根换热管简化为刚性圆柱,并采用双自由度的质量-阻尼-弹簧模型对换热管涡激振动力学模型进行描述,则流向和横向的控制方程分别为:

(3)

(4)

对上述方程采用HHT-α方法[22]进行求解,该方法比传统的Runge-Kutta方法和Newmark-β方法具有更优的算法稳定性和数值耗散性能等。以式(4)为例,进行方程推导:

(1-η)kyt+Δt+ηkyt=(1-η)Fdt+Δt+ηFdt

(5)

(6)

(7)

最终得到关于yt+Δt的方程:

(8)

其中:

(9)

(10)

由式(8)可得:

(11)

将式(6)～(11)编写成C语言程序,以实现单根换热管涡激振动计算中对结构动力学部分的求解。

1.3 几何模型

在固定圆柱绕流计算和单根换热管涡激振动计算中,坐标原点位于换热管中心,流向为x坐标,横向为y坐标,展向为z坐标。二者计算域的相对大小相同,但直径D不同。换热管涡激振动计算的计算域如图1所示,宽为20D,长为30D。基于文献[4,7,9,23]对不同展向长度的圆柱绕流问题进行的研究,选取展向长度为πD。圆柱距离入口10D,距离上下边界10D。入口为均匀速度入口条件,上下边界以及展向的前后边界均为对称边界,出口为压力出口,圆柱表面为无滑移壁面。

图1 涡激振动计算的计算域示意Fig.1 Schematic diagram of the computational domain for vortex-induced vibration calculation

2 网格敏感性分析

表1 敏感性分析的网格模型设置参数

图2 G3网格模型xoy截面示意Fig.2 Schematic diagram of the xoy section of G3 model

为了得到圆柱绕流问题的合理统计数据,根据Franke等[24]的研究,统计时间应大于40个漩涡脱落周期,因此本文的总计算时间约为100个漩涡脱落周期,选取稳定后的70个漩涡脱落周期数据作为有效数据进行统计。

表2 网格敏感性分析结果对比Table 2 Comparison of grid sensitivity analysis results

4套网格的圆柱表面时均壁面速度梯度系数与实验数据[28]的对比如图3所示。壁面速度梯度系数VG定义为VG=|ζ|/Re0.5,|ζ|为无量纲速度梯度的绝对值。可以看出,4套网格的时均壁面速度梯度系数变化均呈现先增后减的趋势,并在θ约为50°处取得极大值,其中分离点为曲线与横轴交点,随着网格数量的增加,分离点逐渐收敛,曲线也逐渐收敛于实验数据,参数θsep逐渐收敛。但需要注意的是,在θ为180°处,计算得到的时均壁面速度梯度系数约为0,明显小于实验数据[28]。由于圆柱绕流尾涡交替脱落,θ为180°处的瞬时壁面速度梯度系数应呈现周期性变化,作时均处理后其值应近似为0,验证了计算的正确性以及数据统计处理的合理性。

图3 不同网格模型时均壁面速度梯度系数分布Fig.3 Time-averaged wall velocity gradient coefficient distributions for different grid models

4套网格的圆柱尾流区不同截面时均流向速度分布与相关实验数据[3,5,29]及数值结果[5]的对比如图4所示。为了将不同截面时均流向速度分布进行对比,对各个截面的数据进行了偏移处理。可以看出,时均流向速度分布由近尾流区(x/D<3)的“U”型向远尾流区(3

图4 不同网格模型圆柱尾流区不同截面时均流向速度分布Fig.4 Time-averaged streamwise velocity distributions of the centerline beyond cylinder for different grid models

图5给出了G3网格模型的圆柱绕流瞬时流场涡量,图中的涡量等值面采用Q准则[30]表示,用无量纲流向速度u/U0进行染色。可以看出,Re=3 900的圆柱绕流问题具有丰富的流动结构,显示出复杂的漩涡脱落现象。同时,附着在圆柱壁面的剪切层呈现二维结构,而分离区和尾流区呈现出高度三维化的涡体结构,具有明显的湍流特征。

图5 圆柱绕流瞬时三维涡量等值面图(Q=10)Fig.5 Isosurface diagram of instantaneous 3D vorticity of flow around a cylinder (Q=10)

综上所述,对圆柱绕流问题而言,近尾流区流动情况复杂,特别是回流区对网格的要求较高,粗网格会导致剪切层较早分离,回流区长度变短,速度型提前由“U”型转变为“V”型;而远尾流区对网格要求较低,不同网格的速度型与实验值吻合较好,仅在峰值处存在误差。也就是说,近尾流区流场对网格更敏感,因此需要对近尾流区的网格进行加密,才能准确地模拟圆柱绕流现象。因而,认为G3可以准确地模拟圆柱绕流现象,并且满足网格无关性要求,则后续工作将基于G3网格模型开展。

3 换热管涡激振动计算

3.1 结构动力学方程求解方法对比研究

表3 不同求解方法的对比Table 3 Comparison of different methods

振幅响应反映的是算法的耗散性能,即算法导致的能量衰减程度。Kane等[31]的研究表明,随着计算时间的延长,Newmark-β类方法的能量衰减要明显小于四阶Runge-Kutta方法。本质上HHT-α方法是由Newmark-β方法衍生而来,具有更优的数值性能。这直接导致方法A3计算得到的振幅响应更接近实验值,而方法A2次之,方法A1最次。

频率响应反映的是算法的弥散性能,即算法导致的周期延长程度。选取了12个振动周期的横向位移y/D的时程变化曲线,起始时刻相位相同,3种计算方法的对比如图6所示。可以看出,随着计算时间的延长,误差逐渐积累导致方法A1和方法A2的振动周期延长,与方法A3相比,方法A1计算的振动周期延长了6.60%,方法A2计算的振动周期延长了23.2%。也就是说,在长时间进行计算时,方法A1和方法A2将会带来较大的计算误差,最终导致计算结果会偏离物理事实。因此,最终选取方法A3(HHT-α方法)对换热管涡激振动计算中结构动力学部分进行求解。

图6 不同计算方法的横向位移y/D的时程变化曲线Fig.6 Time-history curves of transverse displacement y/D for different methods

3.2 单根换热管涡激振动计算

基于建立的换热管流致振动流固耦合计算模型对单根换热管涡激振动问题进行计算,计算参数与J-W的实验参数保持一致。圆柱直径D=0.038 1 m,质量比m*=m/md=2.6,质量阻尼参数(m*+ca)ζ=0.013,水中固有频率fn=0.4 Hz,工质为水,根据实验数据推算出流体密度ρ=995.4 kg/m3,动力粘度μ=0.000 781 Pa·s。时间步长Δt取为0.005 s,保证每个结构固有周期包括大约500个时间步。计算工况为2

考虑到计算中圆柱周围网格需要吸收圆柱边界的运动,采用局部动网格变形技术处理网格变形。对圆柱周围直径为8D范围内采用O形网格进行网格加密,将其内侧直径为3D范围的网格设置为刚体网格区域(A区域),剩余部分的网格设置为变形网格区域(B区域);加密区以外的网格设置为静止网格区域(C区域)。变形前后的网格如图7所示,可以看出,采用局部动网格变形技术变形后的网格具有较好的网格拓扑结构,网格质量较高,不会出现“负体积”现象,可以保证计算的有效进行。同时,由于只有变形网格区域的网格参与变形计算,可以大大节省计算资源。

图7 xoy截面变形前后的网格Fig.7 Grid of the xoy section before and after deformation

3.2.1 振动响应

将不同折算流速下涡激振动响应(横向和流向振幅响应、横向振动频率响应)与J-W[15]的实验数据以及Li[17]和Gsell[18]的计算结果进行对比,如图8所示。可以看出,计算得到的横向振幅和流向振幅与实验数据吻合较好,并且计算结果要优于Li和Gsell的计算结果,计算可以成功捕捉到实验中的初始分支、超级上端分支以及下端分支,各个分支的折算流速范围与实验吻合较好,其中超级上端分支最大横向振幅为1.48D,略小于J-W的实验值1.50D,对应的折算流速Ur为8.35,与实验值更为接近。横向振动频率响应也存在着明显的3个分支,分别与初始分支、超级上端分支和下端分支相对应,且与实验数据吻合较好。3个分支间以“跳跃”的方式过渡,这种过渡方式与实验变化趋势相同。需要说明的是,在某些折算流速下,横向振动可能存在多频现象,这里只选取了横向振动主频。

图8 不同折算流速下振动响应对比Fig.8 Comparison of vibration response at different reduced velocities

当折算流速较小时,流向振幅和横向振幅相对较小,但在2

图9 Ur=3时流向和横向振幅频谱图Fig.9 Spectrogram of streamwise and transverse amplitude at Ur=3

随着折算流速的增大,横向振幅逐渐增大,但流向振幅逐渐减小。当Ur>4时,振动逐渐进入锁定区间,流向振幅和横向振幅逐渐增大,横向振动频率由漩涡脱落频率逐渐趋近于圆柱水中固有频率。在Ur=5时,横向振幅明显增大,振动进入超上端分支。在超上端分支区间5

当Ur>8.65时,横向振幅和流向振幅迅速减小,振动进入下端分支,横向振动频率跳跃至第3频率分支,该频率分支的频率与圆柱在空气中固有频率相近。在整个下端分支区间8.65

3.2.2 运动轨迹

本文对不同折算流速下涡激振动的运动轨迹进行研究,将计算得到的运动轨迹与实验数据进行对比,如图10所示。运动轨迹左为实验值[15],右为本文计算值。对计算得到的运动轨迹进行缩放处理,但整体形状不发生改变。可以看出,计算得到的运动轨迹形状以及流向位移与横向位移的相位差φxy与实验数据吻合较好。

图10 运动轨迹以及流向与横向位移的相位差对比Fig.10 Comparison of motion trajectory and phase difference between streamwise and transverse displacement

当振动处于初始分支和超级上端分支时,运动轨迹为顺时针方向;而当振动进入下端分支时,运动轨迹转变为逆时针方向。在Ur=2时,运动轨迹类似于一条横线,这是因为此时横向位移较小而流向位移较大。随着折算速度的增加,横向振幅和流向振幅均逐渐增大,在Ur=3时,流向振动处于锁定状态,二者振幅相当,并且流向振动频率约为横向振幅的2倍,运动轨迹表现为标准“8”字形。当振动进入超级上端分支区间,横向振幅和流向振幅继续增大,但横向振幅的涨幅明显大于流向振幅,振动轨迹形状逐渐由瘦长“8”字形向倾斜“8”字形转变。当Ur=7.5时,流向振幅达到最大值,随着折算流速的继续增加,流向振幅开始逐渐减小,而横向振幅却继续增大,振动轨迹形状持续变化。当Ur=8.35时,运动轨迹转变为新月形。当振动进入下端分支时,横向振幅和流向振幅均骤减,但横向振幅要大于流向振幅,运动轨迹又逐渐变为“8”字形,但此时运动轨迹为逆时针方向。

3.2.3 尾迹模式

选取xoy截面(z=0),对涡激振动尾迹模式进行研究,图11给出了不同折算流速下的尾迹模式,其中Ωz为z方向涡量,定义为Ωz=∂v/∂x-∂u/∂y。

图11 不同折算速度下尾迹模式对比Fig.11 Comparison of wake mode at different reduced velocities

可以看出,涡激振动尾迹模式分为“SS”模式、“2S”模式、“2T”模式和“2P”模式。当Ur=2时,由于横向振幅较小,而流向振动占主导地位,振动处于“SS”分支,每个振动周期内圆柱后侧脱落一对对称的反向漩涡,与Jauvtis等[14]进行的可视化实验结果相同。随着折算流速的增加,尾迹模式发生转变,当Ur=3时,每个振动周期内圆柱两侧交替脱落方向相反的2个漩涡,该泄涡模式被称为“2S”模式。

当折算流速继续增大时,振动逐渐进入超级上端分支,在这过渡区间上,尾迹模式由“2S”模式向一种新的泄涡模式“2T”模式转变,即每半个振动周期内圆柱一侧脱落3个漩涡,其中一个漩涡的强度远小于另外2个漩涡的强度,且2个强度较大的漩涡方向相同,如图11(c)所示。当振动完全处于超级上端分支时,尾迹模式完全转变为“2T”模式,如图11(d)和(e)所示。当振动进入下端分支时,尾迹模式发生转变,在每半个振动周期内圆柱一侧脱落一对方向相反强度相近的漩涡,称为“2P”模式。

4 结论

1)在进行固定圆柱绕流问题计算时,圆柱近尾流区(x/D<3)流动参数对网格较为敏感,粗网格会使剪切层提前转捩,导致回流区长度变短,而圆柱远尾流区(3

2)与传统的Runge-Kutta方法和Newmark-β方法相比,本文采用的HHT-α方法在求解结构动力学方程上更具优势,体现在较低的数值耗散、较小的数值弥散以及较高的数值精度,使得计算得到的振幅和振动周期更贴近实验数据。

3)在涡激振动计算中,与其他学者的结果相比,本文建立的涡激振动计算模型可以准确预测双自由度圆柱涡激振动行为,计算可以成功捕捉到初始分支、超级上端分支和下端分支,以及相应的“SS”、“2S”、“2T”以及“2P”尾迹模式,超级上端分支范围和最大横向振幅以及不同分支间的过渡方式与实验数据一致。

由于换热管主要以管束的形式存在于蒸汽发生器中,在今后的工作中,可以将本文建立的精细化换热管涡激振动流固耦合模型应用于管束结构流致振动计算中,以支持蒸汽发生器管束结构的设计与优化工作等。