输电导线结冰数值模拟与影响因素分析

高 鹏, 臧五岳, 李丹煜, 徐 枫, 段忠东, 欧进萍

(1. 哈尔滨工业大学(深圳) 土木与环境工程学院, 广东 深圳 518055; 2. 中国电力科学研究院有限公司 输变电工程研究所,北京 100055)

电力工业作为一项经济发展和人民群众生活的支柱产业,其重要性不言而喻。对于电力传输来说,输电线路是远程输送的主要载体。在实际运行中,输电线路通常会受到环境因素的影响,而导线覆冰就曾给世界各地的输电线路安全运行造成严重影响。如2008年南方地区历史罕见的低温雨雪冰冻灾害天气,造成南方电网区域4 216条输电线路被破坏,直接经济损失达1 100多亿元,灾后重建和修复改造需要资金高达390亿元,受灾人口达1亿多人,对我国国民经济造成了巨大损失[1]。

输电导线覆冰通常是由于空气中的过冷却水滴撞击在导线上冻结[2]。国内外学者对导线结冰展开了一系列研究。早期研究人员试图通过对试验数据以及自然条件的分析和总结,建立一个可以表征气象参数和结冰量关系的公式,Langmuir等[3-5]都做过尝试,这些经验模型考虑的因素都比较单一,尽管从结冰概念上是正确的,但是考虑的因素过少,并不能准确地预测导线结冰。Beard等[6]通过试验拟合得到了水滴运动受到的阻力系数。Langmuir等通过试验将水滴中值体积直径离散为七种不同含量的直径,在数值模拟中应用较多,模拟结果更符合水滴碰撞特性。Fu[7]采用Beard等提出的阻力系数以及Langmuir等的尺寸分布模拟得到与试验吻合的结冰冰形。Messinger[8]提出的结冰热力学模型是目前研究广泛采用的结冰模型。Macklin[9]利用风洞试验测定旋转圆柱体的结冰质量,定义了结冰密度计算过程中的参数。Jones[10]基于该参数提出了结冰密度的计算公式。

国内对于导线覆冰数值模型的研究起步较晚。刘和云等[11]在分析了冻雨降水产生结冰的物理过程后, 提出了一个相对简单的结冰预测模型。孙才新等[12]在热力学模型中考虑了电流焦耳热的作用,进一步完善了结冰模型。郭昊等[13]基于雾凇覆冰的数值模拟方法推导了水滴碰撞特性的估算公式,提出了输电线雾凇覆冰的工程估算方法。梁曦东等[14]提出了输电线路时变仿真模型,可以得到时变的结冰外形和结冰质量。刘春城等[15]给出了导线均匀结冰和非均匀结冰的结冰质量和结冰厚度预测模型。Huang等[16]研究了绝缘子的结冰并针对几何参数提出建议,减少绝缘子上的结冰。清华大学开发的三维结冰软件AERO-ICE[17]利用欧拉方法计算水滴场,以自主开发的修正Messinger模型作为结冰热力学分析模型。沈国辉等[18]采用有限元数值模拟方法研究导线的覆冰脱落问题,进行建模参数的敏感性分析和导线脱冰的参数分析。王昕等[19]利用覆冰导线刚性节段模型外加弹簧及配重模拟覆冰输电线路竖向及扭转运动动力特性,制作了D形和新月形覆冰输电线路的气弹模型。以上这些模型在不断发展的过程中对推动导线结冰预测、覆冰导线舞动及冰灾防治研究工作做出了积极作用。

导线结冰的研究主要有试验和数值模拟两种方法。风洞试验的代价较高,在研究不同工况下结冰时比较耗时耗力。本文采用计算流体动力学 (computational fluid dynamics, CFD)数值模拟方法,基于Fluent软件和动网格技术实现导线结冰冰形的数值模拟,研究不同条件参数对水滴运动和导线结冰特性的影响规律,着重分析了导线振动对水滴的碰撞特性和导线结冰外形的影响。

1 输电导线结冰数值模拟方法

1.1 水滴收集系数计算

本文采用拉格朗日法计算水滴运动轨迹。通过Fluent软件中的离散相模型(discrete phase model,DPM)[20]将水滴看成离散相,对每一个水滴的运动状态进行跟踪,建立水滴的运动方程。

水滴与空气之间存在相对运动,当水滴速度小于空气速度时,空气对水滴的曳力成为驱动水滴运动的主要外力;当水滴速度大于空气速度时,空气对水滴的运动会产生阻碍作用。水滴的运动方程可以表示为

(1)

式中:mw为水滴质量;xw为水滴运动位移矢量;t为水滴运动时间;F1为水滴运动受到的阻力矢量;F2为水滴运动受到的重力矢量;ρa为空气密度;Aw为水滴迎风面积;Cd为水滴阻力系数;ua为空气速度矢量;uw为水滴速度矢量。阻力系数Cd根据Beard等根据球体阻力试验数据,得到拟合函数为

(2)

式中,Re为空气与水滴之间的相对雷诺数,可以采用如式(3)计算

(3)

式中:dw为水滴直径;μa为空气黏度。

计算过程中为了使DPM模型中的阻力系数与式(2)保持一致,采用自定义函数(user-defined function,UDF)修改水滴阻力加载到DPM模型中。

撞击在导线表面的水滴数并不均匀,需要通过局部水滴收集系数来反映导线表面水滴碰撞的分布情况。局部收集系数量化了空气来流中有多少水滴撞击到导线表面,可以代表导线表面的水滴分布情况,是研究结冰的关键参数。对于水滴收集系数的计算采用粒子统计法。如图1所示,A0为水滴释放平面,假设其上水滴均匀分布,水滴个数为n0;A1为A0的一部分,假设其上水滴个数为n;A2为导线表面水滴撞击网格,假设其上的水滴全部来自A1。

图1 粒子统计法示意图

水滴收集系数可以定义为释放平面上水滴所围成的面积与撞击表面水滴围成的面积之比,故可表示为

(4)

由于A0中的水滴是均匀分布的,故

(5)

代入式(5),可得

(6)

式(6)中,水滴释放的位置和水滴个数n0由人为预先指定;撞击网格的面积SA2通过UDF中宏函数F_AREA计算得到,撞击网格的水滴个数n也可以通过UDF获取并通过宏函数F_UDMI统计与存储,该方法同样适用于二维导线模型。

由于空气中水滴尺寸并不是单一的,而是有一定的分散性。为了计算方便,常采用水滴中值体积直径(mean volume diameter, MVD)。为了比较单尺寸和多尺寸的水滴碰撞特性,将水滴按照Langmuir D分布离散得到七种尺寸水滴。Langmuir D的分布如表1所示。

表1 水滴多尺寸分布表

计算多尺寸水滴收集系数方法:首先计算不同尺寸水滴各自的收集系数,该过程与计算单尺寸水滴收集系数的方法相同,即前面所述方法;其次将不同尺寸的水滴收集系数按照体积分数进行加权求和,最终得到多尺寸水滴收集系数,如式(7)所示

(7)

式中:ni为第i种尺寸水滴的体积分数;βi为第i种水滴收集系数;β为多尺寸水滴的收集系数。

1.2 导线表面结冰计算

目前,广泛采用的结冰模型是基于Messinger思想建立的结冰热力学模型。对于任意结冰控制体单元,需要满足质量和能量守恒方程

Mimp+Min=Meva+Mout+Mice

(8)

Qimp+Qin+Qair=Qhtc+Qeva+Qclh+Qout+Qice

(9)

式中:Mimp为碰撞到控制体表面的液态水质量;Min为从上一个控制体流入到当前控制体的液态水质量;Meva为控制体表面蒸发的液态水质量;Mout为离开当前控制体的液态水质量;Mice为控制体内生成冰的质量;Qimp水滴撞击到导线表面带来的能量;Qin为从上一个控制体流入当前控制体的液态水带来的能量;Qair为气流气动加热带来的能量;Qhtc为气流与导线表面对流换热带走的能量;Qeva为液态水蒸发带走的能量;Qclh为液态水结冰释放的能量;Qout为离开当前控制体的液态水带走的能量;Qice控制体内冰存储的能量。

在结冰计算过程中假设导线表面的水滴在导线表面流动,不会发生脱落,其流动方式如图2所示,未结冰的水滴会从驻点处平均分为两份沿着导线表面向后流动,图中mw为未结冰的水滴质量。

图2 导线表面液态水滴的流动示意图

控制体内生成冰的质量Mice的计算需要引入冻结系数f。冻结系数的定义为控制体内结冰的质量与流入控制体内所有液态水的质量之比,如式(10)所示

(10)

引入冻结系数后,结冰质量Mice可表示为

Mice=f(Mimp+Min)

(11)

得到控制体结冰质量Mice后,由式(12)计算出控制体表面结冰厚度

(12)

式中:h为控制体表面结冰厚度;ρice为结冰密度,kg/m3;Δt为结冰时间步长。

国内外学者提出过多种结冰密度计算方法,本文采用的是Jones提出的密度计算方法,其表达式为

(13)

2 二维输电导线结冰模型建立

2.1 流体计算域建立

本文选取导线直径D为34.9 mm,采用SSTk-ω湍流模型对导线流场进行非定常计算,采用压力为积的分离式求解器,压力速度耦合采用压力的半隐式方法(SIMPLEC)求解。离散相采用二阶迎风格式,扩散项采用二阶中心差分格式,时间离散采用二阶隐式格式。利用ICEM软件划分网格,计算域尺寸取为30D×20D。导线中心距离入口边界、出口边界分别为10D,20D,距离上、下边界均为10D。导线周围6D×6D为非结构三角形网格加密区,外围区域为结构化四边形网格,网格总数为5.1×104。在导线周围设置30层渐增的边界层,第一层网格高度为5×10-4D,计算得到导线表面第一层网格y+的最大值为0.52。计算域和网格划分如图3所示。计算域左侧为速度入口条件,水平速度等于来流速度;右侧为出口边界采用自由出流条件;计算域的上、下边界采用对称边界条件;对于导线壁面,采用无滑移边界条件。

图3 计算域、网格划分和边界条件

2.2 导线强迫振动模型建立

本文采用强迫振动的形式模拟导线振动,其振动的位移方程为

y=Acos(2πft+φ)

(14)

式中:y为导线横风向位移;A为导线振幅;f为导线振动频率;φ为初始相位。根据导线的振动位移方程可以推导出速度方程

v=2πfAsin(2πft+φ)

(15)

由于初始相位φ不会对导线的频率及振幅产生影响,所以设定导线初始相位0。将以上速度方程通过UDF进行编程和编译,利用CG_Motion宏实现导线的强迫振动,进而在Fluent中构建导线强迫振动模型,并求解振动导线的水滴收集系数。

导线的强迫振动通过Fluent动网格技术的动态层铺法来实现,将计算域分为A1,A2和A3三个分区(见图3(b))。A3分区为随导线一起运动的加密区,A2分区为可变形区域,A1分区为固定区域,不同分区之间通过Interface边界条件进行连接和数据传递。

3 计算结果分析

3.1 结冰模型验证

3.1.1 水滴收集系数

为了验证本文结冰数值模型的正确性,选择了三个已有算例进行验证分析,三个算例的计算条件如表2所示。

表2 水滴收集系数计算工况

采用本文方法计算得到不同算例下的单尺寸MVD和多尺寸水滴收集系数,如图4所示。结果表明,多尺寸和单尺寸水滴收集系数结果整体趋势非常接近,但在水滴撞击的极限处,可以看出考虑了多尺寸分布后,由于存在尺寸大小不一的水滴,水滴碰撞范围更大,与Case1的试验值更加接近,故采用多尺寸分布对水滴轨迹的计算有一定的作用,同时也验证了水滴计算模型的准确性。

图4 不同工况水滴收集系数

3.1.2 结冰冰形验证

选取某个算例对结冰模型进行验证,算例的计算条件包括导线直径为34.9 mm,来流风速为5 m/s,环境温度为-15 ℃,水滴直径为35 μm,空气含水量为1.2 g/m3。

研究共计算了30 min的结冰冰形,分为三个时间段,每个时间段的间隔为10 min。图5给出了导线各时间段下的结冰冰形。从图5可以看出,后一步结冰是在前一步结冰冰形上发展的,这说明结冰模拟是一个多步长的过程,需要不断更新结冰外形进行下一时间步的计算。图6给出了该工况下结冰30 min后的冰形与试验以及模拟结果的对比。从图6可以看出,本部分模型模拟的结冰形状和现有数值模拟结果及试验结果吻合较好。图7所示为不同结冰时间段导线周围时均流线图和速度等值线图。从图7可以看出,不同时间段导线结冰厚度逐渐增加,导线的外形随之发生变化,尾流时均主旋涡的形状也由圆形逐渐拉长,时均主旋涡中心与导线中心的距离也逐渐增大。

图5 导线不同时间段结冰冰形

图6 导线结冰冰形验证

图7 不同结冰时间段导线周围时均流线图和速度等值线图

3.2 不同风速条件下的计算结果分析

风速对导线结冰的影响主要体现在对水滴收集系数的影响。随着风速增大,水滴动能增大,水滴运动轨迹更不容易受到空气的干扰而发生偏移,导致水滴的撞击区域也随之增大,并且水滴收集系数的峰值也增大。图8给出了来流风速分别为2 m/s,5 m/s,8 m/s,12 m/s和15 m/s时导线的水滴收集系数,图9所示为不同来流风速下的结冰冰形,其中水滴中值体积直径为35 μm,空气含水量为1.2 g/m3,环境温度为-15 ℃。

图8 不同风速下水滴收集系数

图9 不同风速下结冰冰形

从图9可以看出,随着来流风速逐渐增大,导线结冰厚度不断增加,同时导线上、下侧结冰的极限位置沿着前缘驻点向后移动,结冰区域不断增大,使得冰形看起来更为饱满。主要原因是:由于风速的增加导致水滴收集系数增大、单位时间内被导线捕获的水滴数增多,结冰的最大厚度也随之增大;另一方面,由于水滴不容易偏转,使得水滴的碰撞极限增大,导线上、下结冰极限也随之增大。

3.3 不同环境温度下的计算结果分析

图10给出了风速为5 m/s,水滴中值体积直径为35 μm,空气含水量为1.2 g/m3,环境温度分别为-21 ℃,-18 ℃,-15 ℃,-12 ℃,-9 ℃和-6 ℃时的导线结冰冰形,分析不同环境温度对导线结冰冰形的影响规律。

图10 不同温度下结冰冰形

从水滴收集系数的定义中可以知道温度对水滴的运动几乎没有影响,温度主要影响结冰时的热力学过程。由图10可以看出,不同温度下的结冰范围是一致的,但随着温度降低,结冰量和结冰厚度不断增加。这是由于温度降低增加了空气的对流散热,也增加了水滴碰撞到导线表面后升高至结冰温度时吸收的热量,导致导线表面散热加快,从而导致单位时间内的结冰增加,并且在温度较高时,导线表面的液态水并不能完全冻结,随着温度降低,冻结的水滴数越来越多,结冰厚度随之增加。当环境温度到达-15 ℃后,结冰厚度变化很小,这是因为此时水滴完全冻结,结冰基本不受温度影响,温度对结冰厚度的影响主要体现在结冰密度上,密度增加导致结冰厚度相对减少。

3.4 不同水滴直径下的计算结果分析

结冰气候条件下,大部分水滴的直径在40 μm以下。图11和图12给出了风速为5 m/s、水滴中值体积直径分别为15 μm,25 μm,35 μm,45 μm、空气含水量为1.2 g/m3、环境温度为-15 ℃时的导线水滴收集系数及结冰冰形。

图11 不同直径水滴收集系数

图12 不同水滴直径结冰冰形

从图11可以看出,随着水滴直径的增大,导线表面的水滴收集系数也随之增大。从数值上来看,最大收集系数并不是随着水滴直径的增大而线性增大,当水滴直径一直增加时,导线表面的水滴收集系数增长会逐渐变得缓慢。在水滴运动过程中,直径较小的水滴会随着空气一起运动而偏移导线,惯性越小受到气流的影响就会越大。相反,较大直径的水滴具有更大的惯性,受到气流的影响较小,会偏离空气流动的方向,从而碰撞到导线上。因此,随着水滴直径的增大,会有更多的水滴保持原有的运动轨迹从而碰撞至导线表面,水滴收集系数会整体变大,水滴碰撞范围也会变大,上、下碰撞极限也会更远。

从图12可以看出,随着水滴直径的增加,导线表面的结冰区域增大,导线表面的结冰厚度也在增加。随着水滴直径的增加,导线表面前缘部分的结冰区域不断变大,整体冰形向两侧发展。这是因为:当水滴直径变大时,水滴运动过程中受气流影响小,更容易被导线捕捉,进而在导线迎风侧会有更大的结冰区域。

3.5 不同空气含水量下的计算结果分析

空气含水量(liquid water content, LWC)是衡量空气中水滴数量的参数,决定了会发生结冰的水量多少,对导线结冰有较大影响。为了研究空气含水量对导线结冰的影响,选取固定导线直径为34.9 mm,风速为5 m/s,环境温度为-15 ℃,水滴直径为35 μm,空气含水量分别为0.5 g/m3,1.2 g/m3,1.5 g/m3。

从水滴收集系数的定义可以知道空气水含量对于水滴碰撞没有影响。计算所得的导线表面结冰冰形结果如图13所示。从图13可以看出,导线表面的结冰量随着空气含水量的增大而增大,结冰厚度也随之增加。当空气中含水量增大时,空气中水滴数量相对变多,单位时间内碰撞导线上的水滴数增加,最终导致导线表面结冰量变多。导线表面结冰范围并没有明显变化,结冰的上、下极限位置没有发生变化。这也是因为在只改变空气含水量的条件下,空气中水滴碰撞特性并没有发生变化,水滴碰撞的极限位置也不会发生改变。

图13 不同LWC下导线结冰冰形

3.6 不同风向角下的计算结果分析

自然界中导线受到的风向并不是固定的,本节为了研究风向对导线结冰的影响,选取了-20°,-10°,0°,10°和20°的风向角进行计算。风向角指的是来流方向和水平方向的夹角,逆时针方向为风向角的正方向,简单示意如图14所示[23]。导线直径为34.9 mm,风速为5 m/s,环境温度为-15 ℃,水滴直径为35 μm,空气含水量分别为1.5 g/m3。

图14 风向角示意图

图15给出了导线在不同风向角下的结冰情况。从图15可以看出,风向角对导线结冰的冰形基本没有影响,与0°风向角下的结冰情况相似,在面对空气来流的前缘驻点处的结冰厚度最大,从驻点沿着导线两侧向后结冰厚度逐渐减小,并且整体的冰形相似。这说明水滴会顺着来流方向在导线前缘驻点处大量聚集,导致该位置处的结冰厚度最大。

3.7 导线振动下计算结果分析

3.7.1 导线振幅对结冰的影响

导线常处于微风振动状态,其幅度通常小于导线直径,选取导线强迫振动振幅分别为0.1D,0.5D和1.0D,导线直径为34.9 mm,水滴中值体积直径为35 μm,风速为5 m/s,导线振动频率为20 Hz,环境温度为-15 ℃,空气含水量为1.2 g/m3。

图16给出了水滴收集系数随导线振动振幅增加的变化。从图16中可以看出,随着导线振幅的增加,导线表面的水滴碰撞极限增大。这是因为当导线沿着横风向发生位移时,原本静止时没有水滴碰撞的位置会因为上、下位移与导线表面逃逸的水滴发生碰撞,从而增大水滴的碰撞极限;也正因为如此,导线上、下表面的水滴数量增多,水滴收集系数大于静止时同位置处的水滴收集系数。在振动频率固定不变时,虽然导线往返运动的次数不变,但是随着振幅增大,导线每次往返运动的范围增加,导致有更多水滴与导线上、下表面发生碰撞,所以随着振幅增加,导线上、下表面处的水滴收集系数增加。驻点处的最大水滴收集系数下降是因为原本会碰撞到导线前缘驻点上的水滴由于导线的竖向运动而导致水滴在导线前缘驻点处发生偏移,从而导致碰撞到导线前缘驻点处的水滴数量变少,水滴收集系数因此而减小。

图16 不同振幅下水滴收集系数

图17所示为导线表面结冰冰形随着振幅的变化规律。从图17可以看出:相比于静止导线,振动导线的结冰极限位置发生了后移,同时导线上、下表面的结冰厚度增加,这是因为导线上、下表面收集了更多的水滴,结冰量也随之增加;而导线前缘驻点处的最大结冰厚度降低,这是因为该位置处捕捉到的水滴数量变少,冻结变小。随着导线振幅的增加,导线上、下表面的结冰厚度由于捕捉的水滴数量增加而增大;前缘驻点处的结冰厚度由于水滴捕捉数量的不断降低而减小。

图17 不同振幅下导线结冰冰形

3.7.2 不同振动频率对结冰的影响

导线发生微风振动的频率范围约为3～100 Hz,故选取导线振动频率分别为10 Hz,20 Hz,30 Hz,40 Hz和50 Hz五组工况进行模拟,导线振幅为1.0D,其余计算条件与前述研究相同。

图18给出了水滴收集系数随导线振动频率的变化。从图18可以看出,随着导线振动频率的增加,导线表面的水滴碰撞极限增大并且收集系数增大。当导线振幅固定时,随着振动频率的增加,导线往返运动的周期变短,因此导线在同样的时间内做往返运动的次数增多,即导线位移的总路程变多,水滴与导线上、下表面接触次数也变多,使得导线上、下表面捕捉的水滴数量增多,水滴收集系数随之增大;与低频率振动相比,频率越高,在相同时间内导线的位移越大,可以捕捉更远处的水滴,因此导线表面的水滴碰撞极限也变大。而前缘驻点处的水滴收集系数则因为捕捉的水滴数量减少而降低。

图18 不同振动频率下水滴收集系数

图19所示为导线表面结冰冰形随着振动频率发生的变化。从图19可以看出,相比于静止导线,振动导线的结冰极限位置发生了后移,同时导线上、下表面的结冰厚度增加,这是因为导线上、下表面收集了更多水滴,结冰量也随之增加;而导线前缘驻点处的最大结冰厚度降低,这是因为该位置处捕捉到的水滴数量变少,冻结变小。随着导线振动频率的增加,导线上、下表面的结冰厚度由于捕捉的水滴数量增加而增大;前缘驻点处的结冰厚度由于水滴捕捉数量的不断降低而减小。

图19 不同振动频率下导线结冰冰形

4 结 论

本文建立基于拉格朗日法的二维输电导线表面结冰冰形预测的CFD数值模拟方法,通过水滴收集系数和结冰冰形与现有试验和数值模拟结果的对比分析,验证了结冰模型的准确性,分析不同条件参数和导线振动对导线结冰特性的影响规律,得到结论如下:

(1)来流风速增大导致水滴动能增加,更多水滴撞击到导线表面,导线表面的水滴收集系数和碰撞极限范围增大。结冰厚度也随着水滴收集系数的增大而增大,结冰沿着导线上、下两侧发展,冰形更为饱满。

(2)水滴直径增大,水滴质量和惯性增大,水滴运动状态更不容易受到空气干扰,导线表面能收集到更多水滴,因此导线表面水滴收集系数和碰撞极限范围增大。结冰沿着导线上、下两侧发展,冰形更为饱满,同时由于结冰密度的增加使得导线结冰厚度的增长并不明显。

(3)随着环境温度降低,导线表面冻结水滴越来越多,结冰由湿增长转为干增长,结冰厚度随之增大。当温度继续下降(-15 ℃),此时水滴已完全冻结,主要影响结冰密度。当空气含水量增加时,相同时间内导线表面收集到的水滴量增多,结冰厚度随之增加。

(4)导线强迫振动振幅越大,导线的水滴收集系数最大值降低,水滴碰撞极限范围增大,结冰厚度最大值降低,结冰逐渐向导线上、下两侧发展。导线强迫振动的频率越大,导线前缘驻点处碰撞的水滴数越少,水滴收集系数最大值降低,导线上、下侧碰撞的水滴数变多,碰撞极限范围越大,结冰逐渐向导线上、下侧发展,冰形更为饱满。