陈崇格,宁小深,李福庚,马嘉辰,李旭晟,胡健
(哈尔滨工程大学 船舶工程学院,黑龙江 哈尔滨 150001)
循环水槽被广泛用于进行船舶工程等领域的水动力试验已有几十年的历史。动力段是循环水槽的重要组成部分,其中的推进叶轮工作时转动主要引起附近区域壁面发生振动。叶轮转动时前方流场形成低压区,后方流场形成高压区,壁面在流体脉动压力的作用下发生变形,而壁面的变形反过来又会影响管内的流场脉动分布,两者相互影响,导致壁面振动的规律复杂。而且壁面振动还会造成循环水槽结构的安全性问题,因此有必要对其振动规律进行分析研究。
国内外已有一些学者进行了关于循环水槽的研究,一方面,在循环水槽中进行了许多试验,包括回转体试验[1-2]、圆柱的涡致振动试验[3]、推进器空泡性能试验[4]以及船模斜航试验[5];另一方面,对循环水槽中流场的流场稳定性[6-8]与湍流度[9]进行了研究。还有Guo等[10]设计了在2个来流方向上效率均较高的适合于循环水槽的叶轮,可以为后续设计大型循环水槽提供参考资料。
循环水槽管道壁面的振动类似于管道振动,在这方面,代路[11]同时考虑三向位移之间的相互耦合作用和边界条件的任意性,采用二维改进傅里叶级数法建立分析模型,实现了对圆柱壳结构振动的精确求解;张允春[12]基于有限元法和机械振动原理,对管道结构的振动机理进行系统的分析,通过改变弯头的角度,对管道系统进行有限元动力响应分析,结果表明管道系统中弯头角度的大小和数量多少,直接影响管道结构振动。且角度越大,振动越大,数量越多,振动也越大,而改变管道尺寸,合理布置支撑,从而改变固有频率和激振力大小,可以有效减小管道结构振动。
直接针对于循环水槽壁面振动规律的研究较少。对于大型减压循环水槽,彭兴宁等[13]用有限元方法计算并分析了其总体和局部自振频率,结果表明其振动特性为扭转加伸缩,且因为推进叶轮叶频低,所以不易引起共振。于昌利[14]联合FLUENT软件和Patran-Nastran软件,不仅模拟了循环水槽内部的流场,还通过计算得到了固有频率,结果表明推进叶轮的一阶叶频范围内可能引起机械共振。
比较国内外关于循环水槽的研究成果可知,已有的研究主要是用有限元方法分析循环水槽的整体及局部模态,并未考虑其内部流场的影响。针对这个问题,本文进行了循环水槽动力段壁面与其内部流场的双向流固耦合计算,分别对时域与频域结果进行分析,对影响循环水槽动力段壁面振动的影响因素进行了分析。
模型对近壁面区域的边界条件有良好的处理能力,但在自由剪切流区域受来流的影响较大,SST模型用k-ω模型处理近壁面区域的流动,用k-ε模型处理自由剪切流区域的流动,综合了k-ω模型在近壁区计算和k-ε模型在远场计算的优点,并在湍流粘性系数的定义中考虑了湍流剪切应力的输送过程,其方程为[15]:
Pk-ρβ*fβ*(ωk-ω0k0)
(1)
σk=0.85F1+(1-F1)
(2)
(3)
β*=0.09F1+0.082 8(1-F1)
(4)
(5)
(6)
φ的定义为:
φ=max(lratioF,1)
(7)
(8)
(9)
(10)
lLES=(0.78F1+0.61(1-F1))Δ
(11)
式中Δ为考虑的网格单元中心和相邻网格单元中心之间的最大距离。
Star-ccm+软件中用拉格朗日方法计算固体结构在空间和时间中的运动。由质量守恒有[15]:
(12)
式中:M、ρ、V和分别为固体结构初始状态时的质量、密度和体积;M0、ρ0、V0分别为固体结构发生变形后的质量、密度和体积。
由于固体结构的体积变化会导致密度发生变化,所以密度为:
(13)
式中:ρ(V,T)为固体结构体积变化过程中的密度;V(T)为固体结构体积变化过程中的体积;T为当前的温度。
在线性各项同性弹性假设中,设Tref为参考温度,固体的应变和密度分别为:
(14)
(15)
式中:K为体积模量;ε为应变;σm为平均应力;α为热膨胀系数。
固体的运动受柯西平衡方程控制,柯西平衡方程表示连续介质的线性动量守恒。固体的动量守恒方程为:
(16)
式中:u为固体位移;b为单位体积的总体积力;σ为柯西应力张量。
对于十分小的应变,即某一网格节点M处的位移,内力是节点位移的线性函数。使用牛顿迭代方法求解控制方程有:
KMNΔuN=rM
(17)
式中:KMN为刚度矩阵;uN为节点的位移;rM为节点M处的残余应力。
由于固体结构的流固耦合计算是一个动态的过程,所以uN满足:
(18)
残余应力为:
(19)
用于双向流固耦合计算的推进叶轮如图1所示,其主要参数如表1所示。
表1 叶轮的主要参数Table 1 Main parameters of the vane wheel
图1 叶轮网格Fig.1 Model of vane wheel
计算域的设置如图2所示,原点位于桨盘面中心处,左侧为速度入口,来流速度为4.67 m/s,右侧为压力出口,静止流体域半径为0.88 m,旋转域直径为0.78 m。固体区域的材料属性设置如下:杨氏模量为6.8×1010,泊松比为0.34,密度为2 700.0 kg/m3。叶轮转速为600 r/min,叶轮叶梢与壁面间的距离为100 mm,壁面厚度为5 mm。在计算域的旋转域和静止流体域的边界处设置交界面以传递计算数据,通过滑移网格实现叶轮的旋转运动。在壁面区域和静止流体域的边界处设置流体—固体交界面,实现流体区域与固体区域计算数据的传递。流体区域用分离涡k-ωSST模型进行计算,固体区域用有限固体模型进行计算。
图2 计算域设置Fig.2 Settings of computational domain
图3 圆柱坐标系Fig.3 Settings of cylindrical coordinate system
图4 壁面振动监测点Fig.4 Wall vibration monitoring point
为了保证数值模拟方法的可靠性,进行网格收敛性分析来保证模拟结果的网格无关性。计算域不同尺度的网格划分详细参数如表2所示。
表2 计算域的网格设置Table 2 Grid settings of computational domain
计算收敛后叶轮推力和转矩的平均值如表3所示,壁面区域监测点的径向挠度变化如图5所示,平均值如表4所示。中尺度网格与细网格相比,推力平均值的相对误差为0.62%,转矩平均值的相对误差为0.86%,径向挠度平均值的相对误差为3.39%,所以可以认为数值模型是准确的。
表3 叶轮的水动力Table 3 Hydrodynamics of the vane wheel
表4 监测点的平均挠度Table 4 Average deflection of monitoring points
图5 监测点径向挠度曲线Fig.5 Radial deflection curve of monitoring point
为了能更清晰地反映壁面的振动规律,将叶轮叶梢与壁面的间隙(以下简写为间隙)分别设置为20、60和100 mm,壁面厚度设置为2 mm,其他参数不变,叶轮转速仍为600 r/min,进行计算和分析。为了便于分析,由于叶轮的转速为600 r/min,故定义叶轮的轴频为f= 10 Hz,叶轮叶数为9叶,则叶轮的叶频为F= 90 Hz。
间隙不同时桨盘面处叶轮正上方壁面区域监测点的径向挠度时域曲线如图6所示,该时间段内径向挠度的平均值如表5所示。间隙从100 mm减小到60 mm时,径向挠度及其波动变化缓慢,但间隙从60 mm减小到20 mm时径向挠度及其波动迅速增加。
表5 间隙不同时壁面监测点的径向挠度Table 5 Radial deflection of the monitoring point on the wall with different gap
图6 间隙不同时壁面区域监测点的径向挠度Fig.6 Radial deflection of the monitoring point of the wall area with different gap
经计算,在间隙为60 mm和100 mm时,挠度时域曲线的小周期在0.011 s,等于叶频对应的周期,但基本上不能看出振动的大周期;而当间隙为20 mm时,能够很明显地看到振动的大周期和小周期,周期虽然有变化,但振动的大周期基本都在0.133 s左右,小周期在0.012 1 s左右。因为当间隙较大时,有足够多的流体缓冲叶轮通过流体施加给壁面的径向力;而当间隙很小时,没有足够多的流体缓冲叶轮通过流体施加给壁面的径向力,壁面附近流场压力较大,导致径向挠度及其波动都较大。
间隙不同时桨盘面处叶轮正上方处壁面区域监测点径向挠度的频域曲线如图7所示,当间隙从20 mm增加到60 mm时,径向挠度波动峰值迅速下降,但从60 mm增加到100 mm时,挠度波动峰值就下降的十分缓慢了。因为壁面固体结构区域的材料属性(密度、杨氏模量和泊松比等)和形状均未发生变化,只是随着间隙的增加,流体与壁面的摩擦力及壁面附近的流场压力都减小了,所以波动峰值下降,但3种情况下都在40 Hz处达到一个峰值,40 Hz很有可能对应于此时壁面的固有频率。
图7 间隙不同时壁面区域监测点的径向挠度波动Fig.7 Radial deflection fluctuations of the monitoring points in the wall area with the different gap
当间隙为20 mm时,径向挠度波动曲线最明显的3个峰值位于7.5、80和90 Hz处,分别对应于一个略小于轴频10 Hz的频率、略小于叶频90 Hz的频率和叶频90 Hz,说明间隙足够小时,由于叶梢和壁面之间流场在径向没有足够的距离缓冲叶轮旋转和壁面振动对流场的影响,振动规律比较复杂,壁面的径向振动还存在一个新的略大于叶轮旋转周期0.1 s的大周期和略大于叶轮叶片旋转周期0.011 s的小周期。但叶频90 Hz对应的挠度波动的值仍然是一个峰值。
当间隙足够小时,径向挠度波动主要由轴频对应的峰值决定,即壁面的振动主要是叶轮旋转大周期内的振动,而当间隙逐渐增加时,径向挠度波动主要由叶频对应的峰值决定,即此时壁面的振动主要是叶轮每个叶片旋转小周期内的振动。
壁面的振动非常微小,且理论上壁面厚度越厚,振动就越小,为了尽可能地避免数值计算误差带来的影响,本节将间隙设置为为20 mm,壁面的厚度分别设置为2、5和8 mm,叶轮转速仍为600 r/min,其他参数不变,进行计算和分析。
循环水槽壁面厚度不同时桨盘面处叶轮正上方壁面区域监测点的径向挠度时域曲线如图8所示,该时间段内的平均值如表6所示。
表6 壁面厚度不同时壁面监测点的径向挠度Table 6 Radial deflection of the monitoring point of the wall with different wall thickness
图8 壁面厚度不同时壁面区域监测点的径向挠度Fig.8 Radial deflection of the monitoring point of the wall area with the different wall thickness
随着壁面厚度从2 mm变化到5 mm时,径向挠度迅速减小,且波动迅速降低;而壁面厚度从5 mm变化到8 mm时,径向挠度的平均值和波动变化趋势减缓,只降低了少量。说明壁面厚度只在较小的范围内增加时,径向的挠度会迅速减小,而当厚度增加到一定的值后,再增大厚度对径向的挠度变化影响很小。
壁面厚度不同时径向挠度的频域曲线如图9所示,当壁面厚度为2 mm和5 mm时,最大峰值均出现在7.5 Hz处,略小于轴频所对应的10 Hz,第2峰值均出现在80 Hz处,略小于叶频所对应的90 Hz,叶频90 Hz处也存在一个峰值。当壁面厚度从2 mm增大到5 mm时,径向挠度波动峰值均迅速降低,而当壁面的厚度再增大到8 mm时,径向的挠度波动均下降得很缓慢,且最大峰值所对应的频率减小到了5 Hz,等于轴频的一半。
图9 壁面厚度不同时壁面区域径向挠度波动Fig.9 The radial deflection of the wall area fluctuates with the different wall thickness
此时壁面固体结构区域的材料属性未发生变化,但随着壁面厚度增加,壁面区域的质量逐渐增大,但壁面附近的流场压力基本不变,因此径向挠度波动减小,波动峰值逐渐下降,但仍在频率为40 Hz处达到一个峰值。
根据3.1节的计算结果,当间隙很小时,较难清晰地反映壁面不受间隙干扰时原本的振动规律。因此本节将间隙设置为100 mm,由于间隙已经足够大,为了能清晰地观测到振动规律,壁面厚度设置为2 mm,叶轮转速分别设置为500、550 和600 r/min,其他参数不变,进行计算和分析。
叶轮转速不同时桨盘面处叶轮正上方壁面区域监测点的径向挠度时域曲线如图10 所示,径向挠度在该时间段的平均值如表7所示。随着叶轮转速的增加,壁面径向挠度逐增加。
表7 叶轮转速不同时壁面区域监测点的径向挠度Table 7 Radial deflection of the monitoring point of the wall area with different rotaional speed of the vane wheel
图10 叶轮转速不同时壁面区域监测点的径向挠度Fig.10 Radial deflection of the monitoring point of the wall area with different rotational speed of vane wheel
叶轮转速不同时桨盘面处叶轮正上方壁面区域监测点径向挠度频域曲线如图11所示,当叶轮转速逐渐增大时,3个方向的挠度波动峰值逐渐增加,但都增加得很缓慢,因为壁面固体结构区域的材料属性和形状均未发生变化,只是随着叶轮转速的增加,壁面附近的流场压力逐渐增加。叶轮转速分别为500、550和600 r/min时所对应的叶频分别为75、82.5和90 Hz,均对应于频域曲线中的一个峰值。
图11 叶轮转速不同时壁面区域监测点的径向挠度波动Fig.11 Radial deflection fluctuations of monitoring points on the wall area with different ratational speed of vane wheel
本节计算时将间隙设置为100 mm,壁面厚度为20 mm,叶轮转速仍为600 r/min,通过改变叶轮模型将叶轮叶数设置为7叶和9叶,和其他参数不变,进行计算和分析。
叶轮叶数不同时桨盘面处叶轮正上方壁面区域监测点的径向挠度时域曲线如图12所示,径向挠度在该时间段的平均值如表8所示。7叶叶轮时壁面径向振动基本上都是叶片旋转小周期内的振动,且振幅较大,而9叶叶轮时壁面径向振动还表现出叶轮旋转大周期内的振动,但振幅小。
表8 叶轮叶数不同时壁面区域监测点的径向挠度Table 8 Radial deflection of monitoring points on the wall area with different number of blades of the vane wheel
图12 叶轮叶数不同时壁面区域监测点的径向挠度Fig.12 Radial deflection of monitoring points on the wall area with different number of blades of vane wheel
7叶叶轮和9叶叶轮桨盘面处叶轮正上方壁面区域监测点径向挠度频域曲线如图13 所示,9叶叶轮叶频90 Hz对应的径向挠度波动峰值远小于比7叶叶轮叶频70 Hz对应的峰值,表明增加叶轮叶数可以有效减低壁面的振动。
图13 叶轮叶数不同时壁面区域监测点的径向挠度波动Fig.13 Radial deflection fluctuations of monitoring points in the wall area with different number of blades of vane wheel
综合本文所有计算过的不同参数的壁面的振动性能,对比图7、图11和图13可知,无论是间隙改变导致流场对壁面的压力变化还是转速和叶数改变引起的叶频变化,虽然在频域曲线上叶频与轴频对应的频率和峰值会发生变化,且有时不明显,但可以发现始终都在频率40 Hz处存在着一个峰值。
结构的固有频率表示固体结构在受到外界激励产生运动时发生自然振动的特定频率,固有频率与外界激励没有关系,是结构的一种固有属性,主要受结构形状和材质的影响。因为壁面区域结构的形状和材料属性均未发生变化,所以一直不变的频率40 Hz很可能是壁面的固有频率。但同时比较可以得出,对于壁面主要的沿径向方向的振动,相比于叶频对应的波动峰值,固有频率对应的峰值要小的多,表明最终影响壁面振动的主要因素还是叶频。
论文结合分离涡模型和有限固体模型,建立了循环水槽动力段壁面、流场与叶轮之间的双向流耦合数值计算模型,分析了间隙、壁面厚度、叶轮转速和叶数对壁面振动规律的影响,研究结果表明:
1) 当间隙足够小时,循环水槽动力段壁面径向振动主要由轴频决定,即壁面的径向振动主要受叶轮旋转大周期内的影响,而当间隙较大时,径向振动主要由叶频决定,即此时壁面径向的振动主要受叶轮叶片旋转小周期的影响。
2)壁面厚度不影响挠度波动的规律,只影响振动的幅值;叶轮转速主要影响挠度波动峰值所对应的频率,但峰值仍与叶轮的轴频、叶频和壁面固有频率有关;
3)增加叶轮的叶数可以非常显著地降低壁面3个方向的振动;虽然壁面区域的固有频率一般都对应于一个径向挠度波动峰值,但峰值都很小,而叶频对应的波动峰值大得多,叶频对壁面径向振动大小起决定性作用。