船用燃气轮机进气部件结冰规律研究

2023-02-15任永鹏万雷王萌王忠义曲永磊

哈尔滨工程大学学报 2023年2期

任永鹏，万雷，王萌，王忠义，曲永磊

(哈尔滨工程大学动力与能源工程学院，黑龙江哈尔滨 150001)

随着燃气轮机应用领域逐渐增多、使用环境逐渐恶劣，燃气轮机进气部件结冰现象逐渐显现[1-2]。当船舶行驶在环境温度为-18～5 ℃的低温海域时，在海风(一般高于9 m/s)及海浪与船体的相互作用下产生大量海水飞沫，飞沫液滴直径在5～5 000 μm。该工作环境极易导致船用燃气轮机进气部件发生结冰问题。而船用燃气轮机进口导叶、导流罩积冰，将导致进气通道面积减小、叶型型线发生改变，影响气动特性，导致压气机压比和效率下降。严重时，进气部件上的冰层脱落会对压气机叶片造成机械损伤，影响燃气轮机的安全运行[3-4]。

国外学者对于燃气轮机进气端部件结冰机理的研究起步较早。Hamed等[5]提出了一种模拟过冷水滴在航空发动机内运动轨迹的三维方法，利用欧拉-拉格朗日法模拟了流场对颗粒运动轨迹的影响。Veillard[6]提出了一种三维准稳态数值模拟方法，对压气机转子、静子以及期间的相互作用进行模拟，提出了一种结合有限元插值法和周向平均法的单元混合交界面的方法。Norde等[7]采用欧拉法，计算冰晶在静子-转子-静子叶栅中的运动轨迹，计算模拟了冰晶蒸发、升华、融化的过程，并将其部分或全部的附着在叶片表面，比较不同尺寸冰晶的总含水量、融化率和撞击质量通量。国内学者主要集中在翼型的结冰问题上。张丽芬[8]利用欧拉-拉格朗日法计算空气-液滴两相流，结合溢流水的溢流方向以及空气速度矢量，整合出一套计算三维翼型表面积冰的方法。苏长明[9]建立了考虑水膜蒸发的三维结冰数学模型。还有部分研究人员在燃气轮机进气道对液滴分布的影响规律方面展开相关研究。屈靖国[10]利用欧拉-欧拉法对亚声速蛇形进气道进行了空气-过冷水滴两相流场的计算，获得进气道出口的流场分布特征。并对三维发动机进口段的支板和整流帽罩进行耦合积冰计算。李静[11]提出了一种三维进气道内结冰参数分布的分析方法。采用欧拉-拉格朗日法获得进气道内表面的水滴撞击特性以及进气道内水滴运动轨迹。

目前国内外对于翼型、航空发动机叶片的结冰问题展开了大量研究。但对于船用燃气轮机叶片结冰，尤其是针对L型进气道影响下的进气部件结冰规律问题研究相对较少。此外，缺乏部件结冰后，气动特性的影响分析研究，以及燃气轮机进气部件结冰临界温度预测研究。

本文利用商业软件FENSAP-ICE对结冰条件下，进气滤清器后的二维船用燃气轮机进口导叶和导流罩表面的结冰情况进行数值模拟，对比分析环境总温、气流流速、液滴平均容积直径、液态水含量等主要因素对冰形的影响规律。同时对典型积冰的叶型和导流罩进行流场计算，对比分析结冰前后叶型、导流罩周围流场变化情况，总结进气部件积冰的影响规律。同时，结合BP神经网络对船用燃气轮机导流罩表面发生结冰的临界温度值进行预测，为后续开发防冰装置和结冰预警装置提供发展方向。

1 数学模型与方法验证

结冰过程的数值模拟一般分为3个步骤：1)求解结冰部件周围空气流场；2)根据空气流场的计算结果进行液滴流场的计算；3)利用结冰计算模型获得部件表面结冰情况。

1.1 数学模型

1.1.1 欧拉法求解液滴撞击特性

欧拉法求解液滴撞击特性，是将液滴看作连续相，引入液滴容积分数的概念后，通过计算液滴的连续方程和动量方程，得到空间网格节点的液滴容积分数和速度分布，进而得到液滴撞击特性的方法。

液滴连续性方程为：

▽·(ρdαud)=0

(1)

液滴动量方程为：

▽·(ρdαuaud)=ρdαK(ua-ud)+ρdαGd

(2)

式中：α为液滴容积分数；ud为液滴速度矢量；ua为空气速度矢量；Gd为液滴重力；K为空气-液滴动量交换系数。

壁面所收集的液滴质量为：

(3)

1.1.2 Shallow-Water结冰模型

液滴结冰模型主要是通过在结冰表面上控制容积的质量守恒和能量守恒，来阐述液滴在结冰过程中热力特性和控制容积中结冰量的方法。结冰热力学计算是从被研究物体的驻点开始，反复分析相邻控制容积的质量守恒和能量守恒，直至所有控制容积分析完毕后，得到每个控制容积内的结冰量。

控制容积内的质量守恒方程为：

(4)

控制容积内的能量守恒方程为：

(5)

(6)

(7)

(8)

式中：hice为冰层厚度；hwater为水膜厚度；uw,in和uw,out分别为溢流流入和流出控制容积的水流速度；hc为对流换热系数；rc为恢复系数；Tl和ul为控制容积边界层上的空气温度和速度。

利用结冰系数可获得各控制容积的结冰量：

(9)

(10)

1.2 方法验证

1.2.1 模型介绍

以0.533 4 m弦长的NACA0012翼型为模型进行结冰方法验证。计算域网格采用O型结构化网格，第1层网格高度0.003 mm，网格膨胀比1.1。以翼型结冰量为参考，对计算模型进行网格无关性验证，最终选择20万网格量进行计算。

1.2.2 边界条件

远场边界：压力远场，气流攻角4°；

壁面：无滑移壁面，结合经验公式给定表面等效粗糙度初始值[12]；

结冰条件：液滴平均直径MVD=20 μm，液态水含量LWC=1.0 g/m3，气流速度ua=67 m/s，喷雾时间6 min，环境总温分别为267.59 K、247.04 K；

结冰表面粗糙度模型：beading模型。

采用时间多步长法进行结冰计算，以1 min为时间间隔计算获得翼型的结冰冰形[13]。

1.2.3 冰形对比

图1中c为弦长，x、y为翼型的横、纵坐标值，y/c和x/c为无量纲量[12-14]。由图可见，计算冰形和试验冰形的吻合程度较高。计算冰形在环境温度较低时，获得的毛冰冰形与试验冰形型线基本一致，能够准确的描述出毛冰的形状特点。在环境温度较高时，冰形整体吻合程度低于毛冰，但计算冰形仍可以较为准确的勾勒出明冰的外形特点。因此，本文所使用的结冰计算方法可用于后续研究。

图1 不同环境温度翼型结冰计算结果Fig.1 Calculation results of airfoil icing under different ambient temperatures

2 船用燃气轮机进气部件结冰

2.1 进口导叶

2.1.1 模型

虽然压气机叶片流动具有三维性和非定常性，但本计算中，不考虑叶根和叶顶等端区位置典型的三维流动特征，忽略导叶沿叶高方向扭转对结冰形状的影响。选用船用燃气轮机2.5级压气机进口导叶50%叶高位置的叶型作为研究对象，即将计算域聚焦为二维。计算域网格划分采用O型结构化网格，近壁面处第1层网格高度0.001 mm，膨胀比1.1。

针对模型进行网格无关性验证，在MVD=10 μm、LWC=0.25 g/m3、环境温度T∞=273 K、空气流速ua=100 m/s、积冰时间t=60 s的条件下对叶型的结冰量进行计算，最终选择20万网格数进行后续研究。

2.1.2 边界条件

远场：给定压力远场边界，气流攻角0°；

壁面：无滑移边界，结合经验公式给定表面等效粗糙度初始值[12]；

结冰条件：该进口导叶模型实际工作中最大进口流速约为130 m/s，参考目前国内外对于海洋环境水雾、盐雾参数，以及船舶进气滤清器分离效率的研究结果[15-17]。滤清器出口处液滴直径一般为5～30 μm，液态水含量一般为0.3～3.5 mg/m3。以滤清器出口滴谱分布中占比较大的液滴尺寸作为边界参数。通过增大LWC，减少结冰时长，保证仿真计算的结冰条件和实际船舶进气道内的气流、水雾环境具有相近的截水常数(Wi=ua×LWC×t)，实现以较短的结冰时长模拟进口导叶结冰的目的，结冰边界条件如表1所示(时间为60 s)。

表1 叶型结冰边界条件Table 1 Blade profile icing boundary conditions

结冰表面粗糙度模型：beading模型。

2.1.3 结冰情况影响

由图2可见，随着气流速度的增大(从工况1增大到工况2)，叶型前缘形成一处厚度较大的楔形冰。这是因为增大的空气流速提高了叶型表面对流换热量，当液滴撞击到叶型前缘后没有立即冻结，而流向叶型的压力面一侧，逐渐冻结、放热，最终形成了向外凸起的角冰。同时，气流速度的增大导致相同时间内，撞击到叶型表面的液滴数量增大，因此叶型压力面上结冰量有所增加。当MVD增大到15 μm(工况3)时叶型的结冰冰形与MVD=10 μm(工况1)时相似，但在积冰区域上有所增大。这是由液滴的气动阻力随MVD的增大而减小，气流对液滴运动轨迹的影响减小，导致大尺寸液滴的撞击面积增大。当LWC为0.5 g/m3时(工况4)，冰形在叶型前缘形成了明显上翘的角冰，角冰的范围从结冰上极限处一直延伸到叶型压力面侧距前缘约10%弦长处。产生该现象的主要原因有以下2点：1)随着LWC的增大，叶型表面的液滴撞击量也随之增大，导致叶型更大的结冰量、更厚的冰层。2)随着LWC的逐渐增大，固壁表面等效粗糙高度会随之增大，导致叶型表面的对流换热增强[12]。当环境温度为276 K时(工况5)，叶型表面仍会发生结冰现象，且冰形为典型的明冰。向上翘起的角冰与吸力面的夹角基本呈90°关系，但滞止区内的冰层厚度相比低温时，变薄明显。形成以上现象的原因是由于高速气流通过叶型时，会对周围的空气产生一定的热力学温降。但由于环境温度整体偏高，因此撞击到叶型表面的液滴不会全部冻结，剩余的液滴会随气流作用最终脱离叶片。

图2 叶型结冰冰形Fig.2 Icing shape of blade profile

2.1.4 结冰对流场的影响

获得冰形后，对其进行网格重构，对比分析结冰前后叶型气动特性的变化[18-22]。此处主要针对T∞=276 K时，即工况5条件下的典型冰形进行流场分析，所设边界条件与未结冰时相同。

从图3可以看出，结冰后叶型进气气流角减小，叶型前缘向上翘起的角冰，使得叶型前缘两侧压力差减小，引起结冰后叶型的进口气流角减小。此外，叶型结冰前后，在叶型的吸力面侧及尾流处都存在一个低速区，该区域内都存在一个顺时针方向、面积较大的分离涡和一个逆时针方向、面积较小的分离涡。这主要由于叶型吸力面侧气流分离导致的，外界气体需要向分离区进行补充，因此产生了一个面积较大的顺时针分离涡。同时顺指针分离涡的下部需要与叶型下方的主流气体进行掺混，相反方向的气流共同作用产生了这个逆时针方向的分离涡。

图3 工况5叶型结冰前后的马赫云图及流线图Fig.3 Mach nephogram and streamline diagram of blade profile before and after icing in case 5

对比结冰前后叶型的流场可以发现，当环境温度为较高的276 K时，由明显上翘的角冰导致叶型吸力面侧气流分离加剧，分离点起始位置从未结冰时叶型弦长的50%处提前至结冰后叶型弦长的1%处，叶型低速区面积也扩大为原来的近2倍，叶型出口气流不均匀度增大。

图4对比了T∞=276 K时，叶型结冰前后总压损失系数的分布和变换情况。本文对叶型总压损失系数ω的计算公式为[23]：

图4 工况5叶型总压损失系数分布Fig.4 Distribution of total pressure loss coefficient of blade profile in case 5

(11)

式中：PT1为叶型进口，即x=-0.02 m处气流总压；PD1为叶型进口气流动压；PT2为叶型出口，即x=0.06 mm处总压。可见环境温度为276 K时，结冰后总压损失系数明显升高，最大总压损失系数较结冰前升高了约22%。同时，结冰后叶型的总压损失系数在尾流中部存在一个先下降后上升的过程。这主要是由于分离区内，2个方向相反的分离涡在掺混处动能有所提升而引起的。此外，结冰后叶型压力损失系数明显增大的区域主要集中在尾流分离区的上部。这表明结冰后叶型的压力损失大小主要受分离区上部的顺时针分离涡影响。分离涡面积增大，叶型总压损失系数也随之增大，导致叶型整合气流的能力下降。

由图5可以发现，结冰后叶型压力面、吸力面的压力梯度明显减小，结冰后叶片所承受的负荷降低。此外，在结冰叶型压力面前缘处压力系数骤减，这主要是楔形明冰对周围气流产生的影响。结冰叶型前缘表面压力的突变可能会导致积冰的破碎或脱落，若脱落的冰块随气流进入压气机内部，则会对压气机叶片产生机械损伤。

图5 工况5叶型表面压力系数分布Fig.5 Distribution of surface pressure coefficient of blade profile in case 5

2.2 导流罩

2.2.1 模型

本节对带有进气道结构的进口导流罩缩比模型进行二维结冰计算，选取中心对称截面为研究对象，忽略三维进气道横向结构对导流罩进气流场的影响。为保证从导流罩出口流出的气流和液滴能够充分发展，不发生回流现象，对导流罩出口流道进行延长处理，延长的长度约为导流罩长度的10倍。

对上述导流罩计算模型进行结构化网格划分，在导流罩处采用C型结构化网格，在喇叭口处采用L型结构化网格，设置导流罩、喇叭口近壁面处第1层网格高度0.001 mm，网格膨胀比1.1。

针对此模型进行网格无关性验证，在空气流速度为100 m/s、MVD=10 μm、LWC=1.0 g/m3、环境温度为273 K、积冰时间为600 s条件下，对导流罩的结冰量进行计算，最终选择38万网格数进行后续导流罩结冰的数值计算。

2.2.2 边界条件

进口：给定压力进口边界，气流攻角0°；

出口：流量出口边界；

壁面：无滑移边界，结合经验公式给定表面等效粗糙度初始值[12]；

结冰条件：导流罩结冰条件设置原则与导叶相同，但由于导流罩尺寸较大，表面的结冰强度相比导叶要小的多，为了能够获得差异明显的结冰冰形，增长了结冰时间，导流罩结冰边界条件见表2(时间为600 s)。

表2 导流罩结冰边界条件Table 2 Center deflector icing boundary conditions

结冰表面粗糙度模型：beading模型。

2.2.3 结冰情况影响

图6展示了5种工况下导流罩的结冰情况及局部放大图。随着气流速的增大(从工况1到工况2)，导流罩表面的冰层厚度明显增加，同时在导流罩滞止点两侧冻结形成2个凸起的角冰，且冰层与导流罩型线间的夹角明显小于90°。这是由于流速的增快增强了冰层表面的对流换热，液滴在滞止点两侧发生溢流而产生的。当MVD增大到20 μm(工况3)时，导流罩表面结冰区域的面积明显增大，且在积冰上、下极限处形成了凸起的角冰，冰层与导流罩型线夹角呈锐角。随着LWC的增大(从工况1到工况4)，导流罩表面结冰冰层的面积、厚度同样有所增大。在环境温度为276 K时(工况5)，可明显地观察到液滴在导流罩表面的溢流情况，冰层呈现滞止点处薄，结冰上、下极限处厚的冰形。整体的冰层厚度相比低温条件下有所减少。但此时的结冰区域面积却比环境温度为273 K时扩大了约19%。

图6 导流罩结冰冰形Fig.6 Icing shape of center deflector

此外，由于船用燃气轮机导流罩前端连接L型进气竖井，该结构导致导流罩积冰发生在其中线以上的滞止区附近。

2.2.4 结冰对流场的影响

将结冰后导流罩进行网格重构，计算其周围流场，分析结冰对导流罩进气整合特性的影响。本节主要针对气流速度为130 m/s和MVD=20 μm，即工况2、3条件下的典型冰形进行流场分析，此处边界条件与未结冰时设置相同。

从图7、8中可以发现，在未发生结冰时，流道上壁面处存在一个面积不大的分离区，此处的气流分离主要是由于高速气流发生了90°的偏转而产生的。此外，从结冰后流场图可以看出，导流罩表面凸起的冰层导致气流发生分离。当积冰较为严重时，燃气轮机上流道的流通面积有所减小，上壁面处的气流分离区消失，主流上移。

图7 工况2导流罩结冰前后的马赫云图及流线图Fig.7 Mach nephogram and streamline diagram of center deflector before and after icing in case 2

图8 工况3导流罩结冰前后的马赫云图及流线图Fig.8 Mach nephogram and streamline diagram of center deflector before and after icing in case 3

从图9也可以发现，工况2条件下导流罩出口，即x=-0.826 m处，结冰前后的速度分布无明显变化。而工况3条件下，结冰导致导流罩出口处下流道速度整体升高，上流道速度略有下降，且主流最大速度区域沿径向偏移了约0.3 m，导致进入压气机的速度场产生了一定的变化。

图9 工况2和工况3条件下导流罩出口处速度分布Fig.9 Velocity distribution at the outlet of center deflector in case 2 and case 3

3 导流罩结冰临界温度预测

3.1 模型建立

采用最小二乘法的方式获得船用燃气轮机进气部件结冰的临界温度计算值，结合BP神经网络进行临界温度的预测[24-25]。具体流程如图10所示。以结冰量是否等于0为标准判断此时的环境温度是否达到结冰临界温度，并建立数据库。随后，利用BP神经网络对初始数据进行训练。最后，对所需条件的导流罩结冰临界温度进行预测。

图10 结冰临界温度预测流程Fig.10 Flow chart of icing temperature prediction

由于本预测网络需要输入气流速度、MVD、LWC 3个自变量，输出结冰临界温度一个因变量。因此网络输入层和输出层神经元个数分别为3和1，经代码调试，并结合神经元个数经验公式：

(12)

式中：l为神经元个数；m表示输入层神经元个数；n表示输出层神经元个数；a为1～10的常数。

最终确定网络层数3层，隐含层神经元个数11个。输入层和隐含层使用了较为常用的双曲正切型函数，输出层也使用了该层网络较常用的线性激励函数。训练函数选用了带有动量和自适应学习速率反向传播的梯度下降函数。

3.2 边界条件及数据库

结冰临界温度初始数据库的建立采用与前文相同的方法对导流罩表面的结冰量大小进行计算。

3.2.1 边界条件的设定

气流速度的变化范围为70～130 m/s；MVD的变化范围为5～50 μm；LWC的变化范围为0.5～2.0 g/m3。以尽量选取边界点，冰形变化点为原则，确定数据库建立数值模拟边界条件，具体参数如表3所示。

表3 结冰临界条件参数与初始计算误差分析Table 3 Error analysis of icing condition parameters and initial calculation

3.2.2 结冰临界温度数据

各工况点的结冰临界温度值如表3所示，以此建立船用燃气轮机进口导流罩表面结冰临界温度初始数据库。

3.3 结果分析

以初始数据库的7组数据和另加5组新数据对该网络计算精度进行检测。初始数据的计算误差为0，如表3所示，5组新数据的计算误差如表4所示(速度均为100 m/s)。

表4 结冰临界温度预测结果误差分析Table 4 Error analysis of icing temperature prediction results

从表中可见，该结冰临界温度预测网络对原始数据的临界温度预测效果较好，对于5组新的工况点，预测相对误差在0.21%以内。可见，在样本点结冰条件范围内，该结冰临界温度预测网络可对导流罩表面结冰临界温度数值进行初步预测，获得误差在正负0.55 ℃范围内导流罩开始结冰的总温值。可为后续船舶进气系统结冰发生条件研究以及防除冰装置设计研究提供方法指导和基础数据。