基于数值计算方法的水下合成射流推力特性

贾连超 耿令波 胡志强 杨 翊 王 超

1.中国科学院沈阳自动化研究所机器人学国家重点实验室，沈阳,1100162.东北大学机械工程与自动化学院，沈阳,110819 3.中国科学院大学,北京，100049

0 引言

传统的水下机器人推进方式主要依靠螺旋桨与舵的配合，当航速较低时很难有效地调节机器人的姿态，并且长时间低速航行会产生低速颤振现象进而产生噪声，机器人隐蔽性较差。为了实现机器人在低速航行时的姿态控制，一种较为常见的解决方案是设置多部喷水推进机构，这又给机器人带来了阻力大、质量增加、系统复杂等不利因素。未来的水下机器人由于工作任务复杂，需要具备低噪声、高隐蔽性、高机动性、结构简单、体积小、水阻小等特点。传统的操纵系统已经不能满足下一代机器人的上述需求，因此，需要一种新型的推进及操纵系统，以满足水下机器人全航速范围内良好操纵性的要求。

合成射流本质上属于非稳态射流。国外研究人员已经通过实验证实，合成射流的推进效率明显高于稳态射流，并且发现射流中涡的出现可以明显提高推进效率[1-8]。WHITTLESEY等[5]的实验结果显示，相比稳态射流，含有涡的非稳态射流使得推进效率提升40%，且与阻力相关的水动力性能提升70%以上。水下的合成射流受到的关注越来越多[9-16]。KRUEGER[1]通过实验发现，涡的产生使得射流推力增大20%以上；KRIEG等[12-13]仿照水母的运动原理设计了一种活塞型合成射流推进器，并对该合成射流推进器进行了一系列推力实验，实验结果显示，虽然合成射流在一个周期内的质量变化为0，但其动量变化不为0，可以产生净推力[12-13]，证明了合成射流用于水下推进、姿态调节的可行性。

合成射流驱动器具有体积小、质量小、结构简单、效率高等优点，且可以完全集成在机器人本体内部，对机器人流体外型影响很小，是一种非常有潜力的新型水下姿态调节方式。目前关于水下合成射流激励器的研究大多集中在实验方面，数值计算方法方面的研究很少，且系统研究网格尺度、边界层以及湍流模型对推力计算精度的影响未见报道。另外，目前对合成射流推力增强性质的机理认识不够清晰。为此，本文研究了水下合成射流的数值计算方法，系统研究了网格尺度、边界层尺度以及湍流模型对推力及流场的影响。建立了合成射流推力的数学模型，并设计了一种出口可偏转的合成射流激励器，通过对该激励器进行数值实验，并结合本文建立的数学模型，阐释合成射流推力增强性质的内在机理。

1 数值计算方法

THOMAS等[17]通过CFX软件对一种采用非对称型激励函数的射流激励器流场进行了数值仿真，并研究了“伴流”对射流诱导动量的作用。KUMAR等[18]在FLUENT软件中采用大涡模拟，对三维射流激励器进行了数值仿真。

合成射流本质上属于一种湍流现象，对湍流的数值计算一般包括直接模拟、大涡模拟及雷诺时均方法模拟，直接模拟与大涡模拟计算量很大，所以在工程应用中普遍采用雷诺时均方法对湍流进行模拟，可以取得理想的精度，因此，本文也采用雷诺时均方法对湍流进行数值模拟。

按照不同的封闭形式，雷诺时均方程又可分为单方程模型(Spalart-Allmaras)、两方程模型(k-ε、k-ω)以及雷诺应力模型，本文主要对比k-ε模型、k-ω模型以及大涡模拟湍流模型对合成射流数值计算的影响。

本文所研究的合成射流激励器的作动部件为活塞，活塞的运动指定为正弦型：

y(t)=asin(2πft)

式中，y为活塞位置；a为活塞振幅，取5 mm；f为活塞振动频率，取10 Hz。

为了尽可能接近物理实际，本文采用动网格对活塞运动进行模拟，网格更新规律如下：

式中，u为流体的速度向量；ug为动网格的速度向量；Γ为扩散系数；SΦ为Φ的源相；A为控制体的面积向量；ρ为密度。

本文所研究的激励器外形如图1所示，激励器的几何参数见表1。本文采用基于动网格的有限体积法进行数值计算[19]，流场边界条件设置如图2所示，激励器底部设置为动边界，激励器的其他边界均设置为无滑移壁面条件。与激励器出口相对的外流场边界设置为压力出口，其余边界设置为压力入口，入口压力与出口压力为环境参考压力，其值为零。

图1 激励器的结构和参数Fig.1 Structure and parameters of actuator

mm

图2 边界条件Fig.2 Boundary conditions

流场设置为黏性、不可压、非定常，流体设置为水，空间、压力、能量均采用二阶精度离散。为保证求解精度，每个周期划分为200个时间步，每个时间步长迭代50次，收敛准则为残差10-5。

为了保证数值方法的可靠性，本文分别研究了不同网格尺度、不同边界层尺度以及不同湍流模型对激励器推力以及流场的影响。除了湍流模型对比之外，本文其他计算均默认采用SSTk-w湍流模型。

1.1 网格尺度的影响

本文研究了不同网格尺度下合成射流推力及流场的变化，为了在保证精度的同时减少整体网格数量，在进行数值计算时，整个流域采用4种不同的网格尺度进行划分：激励器出口网格尺度s1、出口边界层网格尺度s3、激励器腔体网格尺度s4、外流域网格尺度s2，其定义如图3所示。

图3 不同流域内的网格尺寸Fig.3 Mesh sizes in different flow fields

对3种不同尺度的网格分别进行了数值计算，不同网格尺度下的推力曲线对比如图4所示，流场对比如图5所示，其中，网格尺度定义为s1-s3-s4-s2。推力通过对激励器所有表面的压力应力进行积分得到，由图4及图5可知，不同网格尺度下推力及流场基本不变，因此，为减小计算量，本文后续计算均采用0.4-0.04-0.8-10的网格尺度。

图4 不同网格尺寸下的推力Fig.4 Thrust at different mesh sizes

图5 不同网格尺寸下的涡量云图Fig.5 Vorticity cloud chart in different mesh sizes

1.2 边界层尺度的影响

图6 不同边界层网格尺寸下的Yplus值Fig.6 Yplus values in different boundary layer mesh sizes

图7 不同边界层网格尺寸下的推力Fig.7 Thrust in different boundary layer mesh sizes

图8 不同边界层网格尺寸下的涡量云图Fig.8 Vorticity cloud chart at different boundary layer grid sizes

为了研究边界层尺度对合成射流激励器推力及流场的影响，本文对比了4种不同边界层尺度网格的推力及流场计算结果，4种不同的边界层所对应的第一层网格的高度h1分别为0.01 mm、0.04 mm、0.08 mm、0.20 mm。4种不同边界层网格所对应的Yplus值在一个周期内的变化曲线如图6所示。推力及流场的对比结果如图7、图8所示，可以看出，推力及流场对边界层网格尺度的变化并不敏感。

1.3 不同湍流模型的影响

本文研究了不同湍流模型对合成射流推力及流场的影响，所研究的湍流模型分别是Realizablek-ε模型、标准k-ε模型、SSTk-ω模型、大涡模拟(LES)模型。4种不同湍流模型所对应的推力及流场计算结果如图9、图10所示，可以看出，不同湍流模型所得到的推力结果基本一致，但流场细节有差异，大涡模拟与SSTk-ω模型得到的流场非常接近，而Realizablek-ε及标准k-ε计算得到的流场在涡的结构上与前两者不同。Realizablek-ε及k-ε模型计算得到的涡强度明显低于SSTk-ω及LES模型的涡强度。

图9 不同湍流模型下的推力Fig.9 Thrust in different turbulence models

图10 不同湍流模型下的涡量云图Fig.10 Vorticity cloud chart in different turbulence models

如果只关注激励器的推力性质，可以采用标准k-ε模型，该模型对边界层要求较低，可以减小计算量。如果要研究激励器的流场性质，就需要采用低雷诺数模型，如SSTk-ω或大涡模拟，由于大涡模拟计算量较大，故本文默认采用SSTk-ω模型。

1.4 与实验结果对比

为了验证本文数值计算方法的正确性，将数值计算结果与KRIEG等[13-14]给出的实验结果进行对比。KRIEG等并没有给出激励器的具体几何参数，只是提到活塞行程为6.35 mm时排出41.5 mL体积的水，根据这一信息推算出激励器的活塞直径L=91 mm，再根据激励器的L/D=4，推算出此时激励器的出口直径D=23.6 mm。采用前文的湍流模型以及网格尺度，对该激励器在活塞振动频率为11 Hz条件下的推力进行了数值计算，并对得到的推力进行正则化处理，并将结果与KRIEG等的实验结果进行对比，结果显示两者具有很好的一致性，如图11所示。

图11 数值计算和实验数据的对比Fig.11 Comparison between numerical and experimental data

2 推力模型

2.1 合成射流推力

由关于水环境下合成射流的相关研究可知，与稳态射流相比，合成射流具有较高的推进效率[1-3]。为了对合成射流的推力性质尤其是推力增强性质给出合理的解释，本文建立了合成射流的推力模型。在建立推力模型时，本文将激励器内部流体作为控制体，如图12所示。具体建模过程详见文献[20]。

图12 控制体Fig.12 Control volume

2.2 稳态射流推力

本文将具有相同质量流量的稳态射流的推力作为基准值，通过将合成射流推力与稳态射流推力进行对比，来比较不同条件下下合成射流的推进效率。为此，本文建立了稳态射流作用力的数学模型。具体建模过程详见文献[20]。

假设激励器排出的流体可以用一个以激励器出口为底面的圆柱表示，则圆柱高度

(1)

式中,v为流经激励器出口的流体的平均速度；t为时间;a为活塞的位移;k为出口直径与活塞面积之比；T为激励器振动周期。

简化得到，净动量变化

ΔH=ρπ3D2fl2/16

式中，l为激励器一个行程排出的流体柱高度；D为激励器出口直径；ρ为密度；f为激励器活塞振动频率。

则相同质量流量稳态射流在一个周期内的平均作用力

(2)

KRIEG等[12-13]将式(6)作为合成射流的平均推力计算公式，由本文的分析可知，该推力为相同质量流量下稳态射流的推力，合成射流的平均推力比此值偏高。

3 L/D对推进效率的影响

合成射流推进效率易受各种因素影响，其中一个重要因素就是L/D。量纲一参数L/D可用作合成射流效率的一个量度，当L/D值为4左右时，合成射流的效率最高；当L/D值偏离4时，合成射流的效率降低[4,7-8]。

研究图1所示的合成射流激励器在不同L/D时的推力变化，为了便于比较，将推力量纲一化：

不同L/D下的量纲一推力变化如图13所示，图13中的基准线为f2，通过比较数据点与基准线的高低，即可知合成射流推力与稳态射流推力的相对大小。由图13可以看出，L/D为4.3或5.8时，合成射流的推力位于基准线之上，其他L/D条件下推力均位于基准线之下。这说明，L/D为4.3或5.8时合成射流的效率较高，推力大于稳态射流的推力；而L/D<4或L/D>6时合成射流效率较低，推力小于稳态射流的推力。

图13 不同行程比下的平均推力Fig.13 Average thrust at different stroke ratios

对合成射流的研究均表明，合成射流的推进效率受L/D影响较大，普遍认为，当L/D在4附近时，合成射流推进效率最高[4,7-8]，此时合成射流的相对推力大于稳态射流，这与本文得到的结果一致。随着L/D远离最优值，合成射流的效率降低，其推力小于相同条件下稳态射流的推力，这也与本文得到的结果一致。

KRIEG等[13]也给出了不同L/D条件下的推力变化，但给出的结果中，随着L/D的增大，推力呈现单调递减的趋势，这与本文得到的结果不一致。且当L/D>9时，增大频率对推力没有影响，KRIEG等给出的解释是空化的影响，本文数值计算时没有考虑空化，这可能是本文结果与KRIEG等的结果存在差异的原因。

4 机理阐释

图14 喷嘴可偏转的激励器Fig.14 The actuator with nozzle deflection

本文对图14所示的激励器在偏角为0°及90°条件下的推力进行了数值计算。激励器主腔体的回转轴与y轴方向平行，因而，当出口不偏转时，射流的喷射方向为y轴方向；当出口偏转90°时，射流的喷射方向为x轴方向。两种条件下，激励器x、y轴方向的推力曲线如图15所示。

(a)偏角为0°

(b)偏角为90°图15 喷嘴在偏角为0°和90°时x和y方向上的推力Fig.15 Nozzles at 0 and 90 degrees, thrust in the x and y directions

由图15可以看出，当出口不偏转时，y向推力不对称度较大，而x向推力可以忽略；当出口偏转90°时，y向推力接近正负对称，此时x向推力不对称度较高。曲线不对称程度越高，说明一个周期的平均推力越大，对激励器的推力增强性质贡献越大。

两种条件下x、y方向的平均推力见表2，可以看出，当激励器不偏转时，推力主要集中在y向，此时，推力增强也集中在y向。而当激励器偏转90°时，推力主要集中在x向，因而，推力增强也主要集中在x向。

表2 喷嘴偏角在0°和90°时x和y方向上的平均推力Tab.2 Average thrust under nozzle deflection angle of 0 and 90 degree

|f|引起合成射流推力增强的机理解释如下：当激励器喷出流体时，激励器外部流体处于相对静止状态，此时激励器喷出的射流相对于外部静止流体具有较大的相对速度，因此，喷射流体会受到外部流体的阻滞作用，将此喷射流体柱视为刚性，则其受力如图16所示，可以看出，外部流体对控制体的作用力方向与激励器对控制体的力的方向相反，则有

图16 激励器喷出时的受力图Fig.16 Force of the actuator ejected

当流体流入激励器时，吸入流体速度相对于激励器内部流体具有一定的相对速度，因此，吸入流体会受到激励器内部流体的阻滞作用，将吸入流体柱视为刚性，则其受力如图17所示，可以看出，外部流体对于控制体流体的作用力的方向与激励器对控制体的作用力同向，则有

图17 激励器吸入时的受力图Fig.17 Force of the actuator inhalation

由此可以看出，当激励器喷出流体时，激励器受到的作用力比基准值偏大，而当激励器吸入流体时，激励器受到的作用力在数值上比基准值偏小，因而激励器作用力曲线整体上有一个向上平移的趋势。

5 结论

(1)网格尺度及边界层尺度对推力及流场的计算结果影响较小。

(2)不同湍流模型所得到的推力结果基本一致，研究激励器的推力特性可以采用k-ε模型，研究激励器流场特性可以采用SSTk-ω模型。

(3)通过对比分析激励器不同驱动条件下的推力变化可发现，当L/D接近4时，合成射流的推进效率最高；当L/D远离4时，推进效率开始降低。

(4)本文根据推力的数学模型，并通过分析揭示了合成射流的推力增强性质来源于外部流体的阻滞作用。

[1] KRUEGER P S. The Significance of Vortex Ring Formation and Nozzle Exit Over-pressure to Pulsatile Jet Propulsion[D]. Pasadena：California Institute of Technology,2001.

[2] MOSLEMI A A, KRUEGER P S. Propulsive Efficiency of a Biomorphic Pulsed-jet Underwater Vehicle[J]. Bioinspiration&Biomimetics,2010,5(3):036003-14.

[3] KRUEGER P S, GHARIB M.Thrust Augmentation and Vortex Ring Evolution in a Fully Pulsed Jet[J]. AIAA Journal,2005，43(4):792-801.

[4] WHITTLESEY R W, DABIRI J O. Optimal Vortex Formation in a Self-propelled Vehicle[J] J. Fluid Mech.,2013,737(4)：78-104.

[5] WHITTLESEY R W, DABIRI J O. Vortex-enhanced Propulsion[J]. Journal of Fluid Mechanics,2011,668(4):5-32.

[6] KRUEGER P S, GHARIB M. The Significance of Vortex Ring Formation to the Impulse and Thrust of a Starting Jet[J]. Physics of Fluids,2003,15(5):1271-1281.

[7] GHARIB M, RAMBOD E, SHARIFF K. A Universal Time Scale for Vortex Ring Formation[J]. Journal of Fluid Mechanics,1998,360:121-140.

[8] KRUEGER P S, DABIRI J O, GHARIB M. Vortex Ring Pinchoff in the Presence of Simultaneously Initiated Uniform Background Co-flow[J]. Physics of Fluids,2003,15(7):L49-L52.

[9] MOHSENI K. Pulsatile Vortex Generators for Low-speed Maneuvering of Small Underwater Vehicles[J].Ocean Engineering,2006,33(16):2209-2223.

[10] THOMAS A P, MILANO M, G’SELL M G, et al. Synthetic Jet Propulsion for Small Underwater Vehicles[C]//IEEE International Conference on Robotics and Automation.Barcelona,2005:181-187.

[11] KRIEG M, PITTY A, SALEHI M, et al. Optimal Thrust Characteristics of a Synthetic Jet Actuator for Application in Low Speed Maneuvering of Underwater Vehicles[C]//Proceedings of the 2005 IEEE International Conference on Robotics and Automation. Barcelona,2005:2342-2347.

[12] KRIEG M, MOHSENI K. Dynamic Modeling and Control of Biologically Inspired Vortex Ring Thrusters for Underwater Robot Locomotion[J].IEEE Transactions on Robotics,2010,26(3):542-554.

[13] KRIEG M, MOHSENI K.Thrust Characterization of a Bio-inspired Vortex Ring Thruster for Locomotion of Underwater Robots[J].IEEE Journal of Oceanic Engineering,2008,33,(2):123-132.

[14] KRIEG M, COLEY C, HART C, et al. Synthetic Jet Thrust Optimization for Application in Underwater Vehicles[C]//Proc. 14th Int. Symp. on Unmanned Untethered Submersible Technology (UUST).Durham,2005:1-10.

[15] THOMAS A M P. Exploration into the Feasibility of Underwater Synthetic Jet Propulsion[D]. Pasadena：California Institute of Technology,2007.

[16] WHITTLESEY R W. Dynamics and Scaling of Self-Excited Passive Vortex Generators for Underwater Propulsion[D]. Pasadena：California Institute of Technology,2013.

[17] THOMAS A M, ABRAHAM J P. Numerical Simulation of Circular Synthetic Jets with Asymmetric Forcing Profiles[J].The Open Mechanical Engineering Journal,2010,4(1):1-7.

[18] KUMAR A, PANDA P K, KUMAR V, et al. Combined Experimental and Numerical Study of Synthetic Jet in Quiescent Flow[C]//37th National and 4th International Conference on Fluid Mechanics and Fluid Power. Chennai,2010:1-11.

[19] 张瑞，姜峰，杨晋，等.基于动网格的液压缸双向流固耦合分析[J].中国机械工程,2017,28(2):156-162.

ZHANG Rui, JIANG Feng, YANG Jin，et al. Fluid-structure Coupling Analysis of Hydraulic Cylinder Based on Dynamic Mesh Technology[J]. China Mechanical Engineering,2017,28(2):156-162.

[20] 耿令波, 胡志强, 林扬,等.不对称激励函数对水环境下合成射流激励器推力影响及其机理研究[J].机械工程学报,2017,57(3):13-22.

GENG Lingbo, HU Zhiqiang, LIN Yang, et al. Thrust Characteristics of Synthetic Jet Actuator Underwater with Asymmetric Forcing Profile and Its Working Mechanism [J]. Journal of Mechanical Engineering,2017,57(3):13-22.