全附体潜航器螺旋桨桨叶损伤故障辨识机理

刘宗凯， 刘满红， 郭正阳， 李占江

(1.南京理工大学 自动化学院， 江苏 南京 210094； 2.南京越博动力系统股份有限公司， 江苏 南京 211162；3.南京水利科学研究院， 江苏 南京 210029)

0 引言

潜航器在水下与不明物体发生触碰的事故时有发生，一旦螺旋桨、尾舵等关键附体发生故障，可能会导致潜航器瘫痪或失事，造成无法挽回的经济损失甚至大量人员伤亡。智能故障诊断是提高水下潜航器安全性的重要技术，它不仅可以对早期故障进行预警，避免恶性事故的发生，还可以缓解设备维护中维修人员不足或过度维修问题。因此开展潜航器水下故障的在线探测与智能诊断研究非常必要。

目前常用的故障诊断方法主要是基于各种卷积神经网络的深度学习法或迁移学习法[1-2]，深度学习其实是一种有效的数据特征提取技术，该方法已经被应用到图像处理[3]、语言识别[4]、机械故障诊断[5]等领域。深度学习算法的精度主要依赖于特征提取的准确性以及故障样本的数量和质量，如样本的相关性、样本的分类质量、样本数据均衡度、样本分类均衡度等因素。智能故障诊断需积累大量的有效样本，但实际中能采集到的工况极其有限，而随着大规模数值计算的普及，精确的数值仿真可以用来模拟各种实际运行故障，并基于此来建立丰富的样本模型数据库，以弥补数据库建立过程中实际样本少的不足。

螺旋桨作用效果与多种因素相关，如周围流场流动状况[6-7]、上下游绕流物体形状[8-10]、螺旋桨外形[11]等。螺旋桨推进所产生的伴随流场会对航行器的流体动力学性能产生显著影响，Shuai等[12]采用雷诺平均Navier-Stokes (RANS)算法分析了四驱船在航行过程中桨叶所受到的不均衡载荷分布和受力尺度问题。 Köksal等[13]基于KCD-193螺旋桨推进器模型，采用RANS和分离涡模拟(DES)算法并搭配k-ω湍流模型，分析了不同空泡成型条件下的流场特征，推演得到了螺旋桨表面的气蚀强度分布。

螺旋桨桨叶的畸变或者损坏会产生一系列问题，包括驱动力衰减[14]、振动频率或幅值变化以及流场分布不均衡引起的压力脉动[15-16]等。为了探测螺旋桨的振动和噪声，Huang等[17]、Kingan[18]等在考虑螺旋桨、船体及其他附体流场相互作用的基础上，数值分析了螺旋桨对周围场强的扰动。也有学者[19-20]针对船舶螺旋桨数值研究了网格尺寸、几何模型精细程度以及边界层网格类型对数值计算结果的影响，同时也分析了叶梢涡、尾涡以及螺旋桨周围涡结构的特征。胡光忠等[21]采用滑移网格技术，利用有限体积法求解流动控制方程，对螺旋桨在不同进速下的推力系数、转矩系数以及螺旋桨表面的压力分布进行了数值分析，并获取了均匀流场中螺旋桨在不同故障工况下的侧向力和监测点处压力脉动波形，建立了脉动压力频谱图。最后对比了螺旋桨在正常工作状态和桨叶折断故障状态下的频谱图与时域图，通过提取时域波形、特征频率、伴随频率和侧向力系数等特征参量，实现对螺旋桨桨叶折断故障的监测与诊断。欧礼坚等[22]对处于黏性流场中的导管螺旋桨进行了数值计算，分析了导管螺旋桨上的桨叶在不同位置发生折断时对流场的影响以及导管螺旋桨表面压力分布等。前人对机器学习下人工智能的故障检测进行了研究，说明良好的样本库对结果正确性十分重要。文献[23-24]也研究了螺旋桨推进下潜航器所产生的伴随流场对绕流物体的影响。

本文以全附体潜航器桨叶损伤类型的故障辨识为研究目标，通过有限元分析方法获得潜航器在3种不同外形(无损、相邻位置两片桨叶50%损伤和相对位置两片桨叶50%损伤)螺旋桨定转速推进下的流场演化、场强特征要素变化等，从流体动力学角度分析桨叶外形与场强要素之间的耦合作用机理，以建立关键要素之间的联系。为复现不同桨叶损伤对潜航器流体动力学特征的影响，本文首先建立了搭载七叶静音螺旋桨的全附体潜航器模型；其次介绍了数值分析方法，包括网格生成、边界条件设置、控制方程求解等；最后结合流场可视化处理，获得了潜航器的力矩波动曲线变化与螺旋桨损伤的相互关系。本文研究成果可为下一步建立样本库，实现基于迁移学习的故障智能诊断奠定理论基础。

1 潜航器模型和计算区域网格划分

1.1 潜航器模型

图1所示计算模型为全附体潜航器模型，该潜航器模型由围壳、艇身、水平尾舵、垂直尾舵、螺旋桨构成。潜航器体全长4.36 m，艇身直径0.508 m，围壳高为0.184 m，围壳截面为NACA0012翼形体，其弦长为0.35 m，水平和垂直尾翼端截面也为NACA0012结构，其弦长为0.125 m；螺旋桨采用类E1619的七叶静音螺旋桨设计，螺旋桨直径为0.18 m. 模型各部分表面积分别为螺旋桨0.023 73 m2、艇身5.940 6 m2、水平和垂直尾舵0.237 2 m2、围壳0.161 2 m2，总表面积S为6.362 7 m2.

图1 潜航器模型图Fig.1 Submarine model

1.2 计算区域和计算网格划分

图2所示为整个计算区域的网格图。流场计算区域为长方形，该计算区域的长、宽、高分别为20 m、10 m、10 m. 划分网格时需要对流场变化较为剧烈的区域进行加密。如图2中计算区域两侧表面所示，在模型两侧水平方向的网格做了加密，同时由上表面可以看出在潜航器的上下游区域也对网格进行了加密。包裹艇身的流场区域单独划分出一个矩形空间，并采用了逐层加密的结构化网格。

(1)

式中：Δy为y轴方向距离变化量，这里取y+=10.根据该取值计算出近壁面第1层网格厚度为1.208×10-4m，网格高度增加比为1.2，网格层数为15层，以此建立边界层网格。图3所示为流场区域网格切面图。由图3可以看出，这里采用O形网格，可使包裹艇身的正方形区域网格最密集。

图3 流场计算区域切面处网格图Fig.3 Grid section diagram of the flow field calculation area

图4所示为潜航器壁面网格分布。在围壳、尾舵和艇身周围均采用C形和O形网格进行离散，并在围壳和尾舵等流体变化较为剧烈的区域采用加密网格进行处理。对于螺旋桨扰动区，由于螺旋桨外形特殊且不停转动，不适合用结构化网格，因此在螺旋桨周围单独划分了圆柱形区域，该区域中对桨毂、桨叶以及圆柱流域均采用非结构化四面体网格进行离散。生成网格后该区域的网格节点数为172 679，网格单元数为886 700. 整个流场区域的网格节点数为1 905 855，网格单元数为2 571 642.

图4 潜航器壁面网格图Fig.4 Grid diagrams of submarine surface

2 数值方法和初始条件

2.1 数值方法

数值模拟主要通过求解三维非定常RANS方程组(URANS)和湍流封闭方程组来实现。在螺旋桨扰动的圆柱形区域内采用动网格技术，因此该区域中的控制方程首先要建立在笛卡尔坐标系{xl}(l=1，2，3，表示该坐标系的3个方向)下，通过曲线浸没边界方法将笛卡尔坐标系下的位置和速度矢量{Uj}(j=1，2，3，表示该坐标系的3个方向)通过(2)式和(3)式转换到广义曲线坐标系{ξi}(i=1，2，3，表示该坐标系的3个方向)下[25]。URANS方程可以表示如下：

(2)

(3)

(4)

δij为单位张量，μt为动态涡黏度，

(5)

(6)

k为湍流动能，β*为闭环系数，d为到最近壁面的距离。

计算时为了使URANS方程组闭环，这里添加了剪切应力传输(SST)k-ω模型。则在广义曲线坐标系{ξi}下，添加了k-ω模型的URANS方程组[26]可以表示如下：

(7)

(8)

(9)

F2为混合函数，

(10)

(11)

三维URANS方程的求解，对流项和耗散项分别采用3阶加权无振荡(WENO)格式和三点2阶中心差分格式(SCDS)格式进行离散。在时间尺度上控制方程是通过Crank-Nicolson方法进行的离散，离散后控制方程由全隐式Newton-Krylov方法进行求解[27]。文献[28]中使用了该方法，其正确性也得到了验证。

2.2 初始条件

速度入口边界条件：来流速度U∞=1 m/s，ρ=1 000 kg/m3，沿x轴正向流入，湍流动能密度为中等湍流密度5%，时间步长为0.01 s.

压力出口边界条件：出口边界流体的平均静压为0 Pa，参考压力为0 Pa.

壁面边界条件：艇身、围壳、尾舵、桨毂、桨叶都设置为无滑移壁面边界条件。

如图5所示，桨叶按逆时针顺序分别编号为1号～7号。根据桨叶损伤特征可以分为3种工况：工况1，完整无损伤；工况2，相邻两个(4号和5号)桨叶损伤50%；工况3，相对两个(2号和6号)桨叶损伤50%.

图5 3种工况下桨叶损伤结构图Fig.5 Structural diagrams of blade damages under three work conditions

3 桨叶损伤与力矩扰动的耦合机理

未来在工程应用中需要对潜航器运行状态进行监控，再将采集到的数据通过滤波和信号处理以辨识桨叶损伤类型，一般是通过加装惯性测量元件(如惯性导航元件和加速度计等)测量艇身姿态和方位，或者在艇身外壁面贴装压力传感器来探测表面压力。为更好地研究不同桨叶损伤对潜航器运动和流场特征的影响，主要针对3种不同的桨损类型，依托流体力学的基本理论对桨叶损伤与力矩扰动的耦合机理进行研究。

3.1 工况1：完整桨叶

图6 艇身力矩波动曲线(工况1)Fig.6 Fluctuation curves of submarine moment in workcondition 1

为了更好地探索图6中波形产生的原因，需要分析桨叶与流场的作用机理，图7给出了不同时刻处无损伤螺旋桨(工况1)推进下艇身尾部流线和壁面压力分布图，其中切面1和切面2是在x=2.5 m和x=3.7 m位置处垂直于x轴半径为0.38 m的圆，反映的是艇身近壁区域法向速度梯度变化特征。由图7可见：切面1与艇身相交汇的区域呈现出淡黄色的环形，该区域为近壁面边界层，该环形区域流体流速较小，同时在切面1上部中轴线位置也出现了一个垂直低速带，这是前方围壳阻碍所引起的；切面2中与艇身交界的环形低速区较切面1有所增厚，这是因为该切面位于艇尾的收缩区域，上游艇身的阻碍引起了该区域速度减小。由图7(a)～图7(e)可以看出，流线流经螺旋桨时由于受到无损螺旋桨叶的推动，该区域流体被加速，在螺旋桨叶梢区域流体速度增加最为显著，同时流线在叶梢区域与桨叶分离，形成螺旋形曲线并加速向后运动。由于螺旋桨后侧流体流速的增加，其后侧流体压力降低，并形成了一个负压区。

图7 不同时刻的艇身壁面压力与切面速度云图(工况1)Fig.7 Submarine surface pressure and velocity section nephograms in work condition 1

图8所示为工况1下螺旋桨及下游区域流线侧视图。该侧视图主要被划分为4个区域，分别为A1侧视区域、B1侧视区域、B1侧视流线轮廓上和B1侧视流线轮廓下。对比不同时刻可知，在A1区域由于叶梢和导边的作用，使得在螺旋桨周围产生了包裹螺旋桨的密集流线，且该区域速度幅值较大，并随着时间的推移，该区域黄色流线分布形态并没有产生显著的变化；图8中B1流线轮廓的上边缘和下边缘，随时间的变化分别产生了上下摆动的形态，并在B1侧视区域内不断有尾涡生成。

图8 不同时刻的螺旋桨及下游区域流线演化图(工况1)Fig.8 Streamline evolutions of flow fields around the propeller and downstream area in workcondition 1

图9所示为工况1螺旋桨下游x=4.46 m位置处切面压力云图。此外由于螺旋桨由7个桨叶组成且呈现非轴对称分布，图9中蓝色负压区随着螺旋桨的旋转绕轴心不断滚动，这种滚动频率较螺旋桨转动频率低，周期大约是1.9 s，对应于图6中俯仰力矩系数CMp,y和偏航力矩系数CMp,z包络大幅度波动的频率，也对应图8中流线上下摆动的频率。这种不稳定影响会在螺旋桨周围均衡出现的，使螺旋桨下游区域压力随时间的变化总体保持不变，因此图6中压差力矩系数的平衡线基本位于零位置。

图9 不同时刻的螺旋桨壁面与下游压力切面云图(工况1)Fig.9 Pressure nephograms of propeller and downstream section in work condition 1

3.2 工况2：相邻位置存在两个50%损伤的桨叶

图10所示为相邻两桨叶损伤时(工况2)整个艇身的压差力矩系数随时间演化曲线。由图10可以看出：相对于完整桨叶工况1，俯仰力矩系数CMp,y的平衡线由零位置产生了向上的平移，表明艇身壁面压差力的变化使得其产生了一个抬头的趋势；偏航力矩CMp,z向负方向偏移，表明艇身产生了一个向右侧偏摆的转动趋势，因此可以推测在损伤螺旋桨的推动下，在艇身左弦后侧区域出现一个负压区。这是因为桨叶损伤产生的不稳定扰动涡使艇身尾部特定区域出现了压力较低的区域，该区域流速增加压力减少，而流速的增加使得该区域压力减少边界层的黏性摩阻增大，因此引起了压差力矩系数平衡位置的上下平移。此外，CMp,y和CMp,z二者之间的相位差φ大约是0.75 s较工况1的0.5 s增大了约50%.

图10 潜航器压差力矩波动曲线(工况2)Fig.10 Fluctuation curves of submarine pressure moment in work condition 2

图11所示分别为不同时刻下相邻50%桨叶损伤的螺旋桨(工况2)艇身壁面压力和尾部流线分布图。由图11可以看出，在螺旋桨中心区域下游附近出现一个约是1/2螺旋桨尺度的涡，该涡向下游扩散，且呈现向上偏移的趋势。

图11 不同时刻的艇身壁面压力与切面速度云图(工况2)Fig.11 Submarine surface pressure and velocity section nephograms in work condition 2

图12所示为工况2下螺旋桨及下游区域流线侧视图。由图12可以看出，相对于工况1的A1侧视区域，工况2中的A2侧视区域并不能完全框住黄色高速区域的流线，在A2区域的下游也出现了黄色流线，尺度略小于A2区域的流线。这主要是因为在折损桨叶的作用下部分高速流体产生于折损桨叶的断面区域，并在B2区域产生了小尺度涡，其与外流场的相互作用使尾流B2区域轮廓上下边界产生了上扬，且不随时间显著地变化。

图12 不同时刻的螺旋桨及下游区域流线演化图(工况2)Fig.12 Streamline evolutions of flow fields around the propeller and downstream area in work condition 2

如图13所示，由于两个相邻桨叶的损伤，整个螺旋桨在旋转过程中其尾部的蓝色负压区产生了一个位于第Ⅲ象限(负y轴和负z轴区域，左弦后下方)的偏心。由于这种压力偏心的存在使得其压力矩系数CMp,y和CMp,z的平衡线分别产生了向上和向下的偏移。相对于工况1的情况，该切片云图中蓝色负压区域并没有随着桨叶的转动产生明显的滚转，而一直出现在图中的第Ⅲ象限，其本质是桨叶转速一定的情况下，旋转过程中连续两个桨叶的缺损在同一位置出现的概率较高，该区域折损桨叶的扰动产生了小尺度涡，进而产生了推力偏心。

图13 不同时刻的螺旋桨表面压力与下游切面压力切片云图(工况2)Fig.13 Pressure nephograms of propeller and downstream section in work condition 2

3.3 工况3：相对位置存在两个50%损伤的桨叶

图14所示为相对位置存在两个50%损伤的桨叶条件下(工况3)艇身壁面压差力矩系数随时间变化曲线。由图14可以看出：压差力矩系数CMp,y和CMp,z波形的平衡位置分别位于0.000 3和0.000 2附近，相对于零位置都略微产生了正向偏移，但偏移幅值较工况2时小；CMp,y和CMp,z波形包络的波动周期大约也还保持在2 s左右，同时二者相邻峰值的相位差φ也基本恢复到了0.5 s.因此压力所产生的俯仰和偏航力矩波形峰值的时间差主要取决于桨叶的对称性，可以推测即便是桨叶损伤但只要损伤位置有一定的对称性，也能保持二者的相位与无损时一致。

图14 潜航器压差力矩波动曲线(工况3)Fig.14 Fluctuation curves of submarine pressure moment in work condition 3

图15所示为不同时刻下在两个相对位置同时存在50%桨叶损伤的螺旋桨(工况3)艇身壁面压力和尾部流线分布图。由图15可以看出，小尺度涡从桨毂后部产生并向上方扩散，其在前方艇身绕流的作用下其向上扩展趋势被抑制，并相互融合向下游移动。

图15 不同时刻的艇身壁面压力与速度切面云图(工况3)Fig.15 Submarine surface pressure and velocity section nephographs in work condition 3

图16所示为工况3下螺旋桨及下游区域流线侧视图。由图16可以看出：在A3区域还是分布了速度较大且较为密集的流线，同时在A3下游螺旋桨桨毂附近的区域，也有橙红色的流线生成，不同与工况2中的A2，这个区域的流线分布较为均匀且并未出现工况2中流线的显著倾斜；在B3区域依然有折损桨叶所产生的小尺度涡，其流线轮廓的上下边缘呈现出一个先上扬后下摆的轨迹，这也是图14中力矩系数平衡线偏移零位置较少的原因。

图16 不同时刻的螺旋桨及下游区域流线演化图(工况3)Fig.16 Streamline evolutions of flow fields around the propeller and downstream area in work condition 3

如图17所示，由于两个相对桨叶的损伤，整个螺旋桨在旋转过程中其尾部的负压区也产生了一个向下(负y轴方向)的偏心，但是不同于工况2，这种偏心基本对称于垂直轴线，并随着螺旋桨的旋转产生左右小幅度的摆动，这种对称分布是该状况下压力矩平衡线偏离零位置较少的原因。

图17 不同时刻的螺旋桨壁面与下游压力切面云图(工况3)Fig.17 Pressure nephographs of propeller and downstream section in work condition 3

4 结论

本文以螺旋桨推进下全附体潜航器为研究对象，基于流体动力学基本理论通过求解三维URANS方程和湍流封闭方程获得均匀来流条件下，具有3种不同桨损特征的定转速推进的全附体潜航器所形成的扰动流场与力矩波动的演化过程，并分析了桨损类型对场强特征要素及其空间分布的影响，进而梳理出桨叶损伤特征、流场演化以及艇身动力学特征之间的耦合作用机理。得出主要结论如下：

对于3种工况，其压差俯仰力矩和偏航力矩系数曲线都会存在两种成分的波动：一个是低频大幅度波动，其平衡线位置与桨叶对称性相关；一个是高频小幅度波动，其频率与螺旋桨转动频率一致。

1)工况1时，由于桨叶绕轴中心对称，艇身俯仰和偏航压力矩系数包络的平衡线都出现在零位置，且相邻峰值的相位差较小。此外完整桨叶的叶梢都位于相同半径的圆周上，使得螺旋桨下游B1区域也没有产生显著的小尺度涡。

2)工况2时，俯仰和偏航压力矩系数包络的平衡线分别产生了正向和负向偏移，同时二者的相位差也被拉大，这是因为在相邻两折损桨叶的作用下桨毂下游产生了尺度较小的偏心涡，其与外流场的相互作用使下游流线轨迹向一侧偏移，使推力出现了偏心。

3)工况3时，由于相对桨叶的缺损使其也产生了小尺度的涡，该涡绕轴线有较好的对称性。俯仰和偏航压力矩系数包络的平衡线又基本回归到了0值附近，但较工况1的平衡位置略有上移，相位偏移也与无桨损时类似。

综上所述，桨叶推力的偏心位置决定了压力所引起的俯仰和偏航力矩系数平衡线相对零位置的偏移量，而二者的相位差则反映了桨叶的对称性，下一步结合黏性力矩的变化，可以推测出桨叶损伤类型。