本发明属于水下仿生机器人技术领域，具体涉及一种仿生波动鳍的推进性能数值计算方法。
背景技术：
智能水下机器人是新一代海洋装备发展的关键技术，对提升海洋军备竞争力有着巨大推进作用，已成为各国政府、科学家及企业高度关注的前沿科技领域。在高性能水下装备的研发中，自然界给予人类重要启示，仿生波动鳍水下推进器就是从鱼类高效运动获得灵感的典型例子。目前，国内外针对仿生推进器外观、结构和性能已有大量研究。相比于传统水下航行器，仿生推进器具有推进效率高，隐蔽性强，体积小，价格低廉等优势，其在资源探测，环境监测，水下作业等多方面都显示出巨大潜力。
自然界中，黑魔鬼鱼依靠在鳍面传递行波从而产生推力，通过改变鳍面行波方向能快速反转运动方向，在低速下稳定性强，机动性强，推进效率高等特点。研究者从自然界黑魔鬼鱼高效的运动方式得到启示，基于黑魔鬼鱼运动机理，开展仿生波动鳍水下推进器研究，以克服传统推进器效率低、隐蔽性差、结构复杂等缺点。仿生波动鳍水下推进器具有功能多样、行动高效、结构灵活、性能优越等特点，在海洋探测、情报收集以及两栖作战方面有着巨大应用前景。
水动力性能分析是仿生推进器结构设计和运动控制的理论基础，对提高推进器综合性能起到至关重要的作用。传统性能测试方法主要包括水洞和风洞实验，存在成本高、耗时长等缺陷。近几十年，计算流体动力学(cfd)技术蓬勃发展，以电子计算机为工具，应用各种离散化的数学方法，对流体力学的各类问题进行数值实验、计算机模拟和分析研究，以解决各种实际问题，极大节省人力物力。在仿生波动鳍研究领域，众多研究者结合cfd技术开展了波动鳍运动机理和流场结构演化规律研究，但针对波动鳍推进力的具体影响因素和变化规律尚未得出具体结论，已有研究并不能对仿生波动鳍推进器的结构设计与优化起到指导作用。
技术实现要素：
本发明的目的在于提供一种仿生波动鳍的推进性能数值计算方法，以克服现有技术存在的不足。
为达到上述目的，本发明采用如下技术方案：
一种仿生波动鳍的推进性能数值计算方法，包括以下步骤：
步骤1)，建立待计算仿生波动鳍的仿真计算物理模型；
步骤2)，确定仿真计算物理模型外流场范围以及边界条件；
步骤3)，根据外流场范围以及边界条件建立粗糙网格并进行波动鳍运动初始化，得到粗糙网格仿真计算物理模型；
步骤4)、对粗糙网格仿真计算物理模型进行精细网格划分得到精细网格仿真计算物理模型；
步骤5)，对精细网格仿真计算物理模型进行仿真参数和运动条件设置；
步骤6)、基于步骤5)设置的仿真参数和运动条件参数计算仿生波动鳍推力：选定总计算时间、时间步长和计算迭代次数，采用迭代法进行仿生波动鳍推力计算，得到最终仿生波动鳍推力性能数值。
进一步的，确定外流场区域范围指确定波动鳍周围流场范围，采用尺寸大小确定波动鳍周围流场范围；波动鳍流场边界条件指流场前端半圆面设置为速度入口条件，后端平面设置为压力出口条件，周围圆柱面设置为壁面条件。
进一步的，对仿真计算物理模型进行网格划分，确定鳍面的初始空间位置。
进一步的，根据鳍面运动方程确定鳍面的初始空间位置，鳍面运动方程如下所示：
其中，参数d和s是鳍面上一点位置的空间坐标，θm是波动鳍运动幅值，t是波动鳍运动周期，φ是波动鳍运动相位差，λ是行波波长，无量纲系数ry和rz用于调节运动坐标；取时间t＝0s，根据方程可以确定波动鳍初始空间位置。
进一步的，对步骤3)得到的粗糙网格仿真计算物理模型进行精细网格划分：将粗糙网格仿真计算物理模型划分为两个计算域，靠近仿生波动鳍模型的网格进行加密，网格类型为四面体网格，设置为一个计算域；远离仿生波动鳍模型的流场，网格类型为六面体网格，设置为另一个计算域；两个计算域通过共节点面连接，精细网格总数不少于200万。
进一步的，将步骤4)得到的精细网格仿真计算物理模型导入软件fluentmeshing中，进行数值求解相关参数设置，采用非定常和压力基的湍流求解器求解流场。
进一步的，采用雷诺平均纳维-斯托克斯rans方程求解流场：
其中，
进一步的，涡粘系数模型采用realizablek-εmodel，壁面函数采用yplus不敏感处理；压力速度解耦算法采用piso算法，时间项采用一阶隐式欧拉格式，压力梯度插值采用二阶格式，对流项采用二阶迎风格式；压力项的松弛因子降低到0.15以提高数值计算的收敛性；其中涡粘系数模型的运输方程为：
其中，gk和gb分别是由于流场平均速度和浮力导致的湍动能生成项，ym代表可压缩湍流中脉动扩张对总耗散率的贡献，σk和σε是湍流普朗特数，c2,c1ε和c3ε是常数，sk和sε是用户自定义的源项。
进一步的，采用分离式求解算法，在初值条件下先对流场速度进行求解，并进行压力矫正，使得求解的物理量满足守恒方程；其次对湍动能标量进行数值计算，在流场物理量求解完毕后，根据求解值对推力进行计算：ft＝t(t)-d(t)
其中t(t)和d(t)分别代表压力和摩擦力产生的推进作用和阻碍作用；
t(t)和d(t)分别为：
其中，p是单位压力，τj1是x方向的粘性应力张量分量，nj是方向向量，a是波动鳍总面积。
进一步的，根据步骤6)得到的仿生推进器运动产生的推力对仿生波动鳍推进性能进行分析：通过改变波动鳍运动的幅值与频率，获取波动鳍推力的变化规律。
与现有技术相比，本发明具有以下有益的技术效果：
本发明一种仿生波动鳍的推进性能数值计算方法，根据仿生学原理建立相关物理模型，确定仿生推进器外流场范围以及边界条件类型，采用了两次网格划分技术，提高了计算准确性，建立粗糙网格并进行波动鳍运动初始化，完成初始化后再基于模型划分精细网格，并采用动网格方法进行仿真计算；本发明提出的数值方法基于粗细两套网格，能够更加精确控制推进器初始状态，提高了仿真结果的准确性，得到了一个周期内波动鳍鳍面压力分布规律，流场随波动鳍运动演化规律以及波动鳍推力与波动频率之间变化规律，有效预测了仿生推进器的推进性能，为仿生推进器结构设计优化和运动控制奠定理论基础。
附图说明
图1为本发明的数值计算方法流程示意图。
图2为本发明的波动鳍周围流场示意图。
图3为本发明的波动鳍运动初始化示意图。
图4为本发明的波动鳍精细网格图。
图5为本发明的推力数值计算流程图。
图6为采用本发明数值方法得到的波动鳍压力分布图。
图7为采用本发明数值方法得到的波动鳍流场漩涡分布图。
图8为采用本发明数值方法得到的波动鳍推力曲线。
具体实施方式
下面结合附图对本发明做进一步详细描述：
一种仿生波动鳍的推进性能数值计算方法，包括以下步骤：
步骤1)，建立待计算仿生波动鳍的仿真计算物理模型：
仿真计算物理模型包括鱼体、左侧波动鳍和右侧波动鳍，两侧鱼鳍与鱼体并未连接。鱼体为主要部分，在计算模型中被简化为一个类流线型的刚体；左侧和右侧的鱼鳍为长条形，尺寸为0.4m×0.15m×0.005m，通过模拟黑魔鬼鱼波动鳍运动，从而产生向前推力；仿生波动鳍是指模仿自然界黑魔鬼鱼运动机制而设计的一种水下仿生推进器。
步骤2)，确定仿真计算物理模型外流场范围以及边界条件；
具体为确定外流场区域范围和边界类型，模拟物理模型周围流场环境；
如图2所示，确定外流场区域范围指确定波动鳍周围流场范围，本申请中采用尺寸大小确定波动鳍周围流场范围，外流场尺寸为6l×4l×4l，其中l是仿生推进器总长度，本申请l取0.5m；确定波动鳍流场边界条件指流场前端半圆面设置为速度入口条件，后端平面设置为压力出口条件，周围圆柱面设置为壁面条件；
步骤3)，根据外流场范围以及边界条件建立粗糙网格并进行波动鳍运动初始化；
如图3所示，波动鳍运动初始化是指，由于波动鳍运动是一个周期运动，鳍面初始时刻所处位置将直接影响仿生推进器的推进性能，因此需要根据相关理论确定初始位置；具体方法是，采用软件fluentmeshing对得到的仿真计算物理模型进行网格划分；由于该套网格不参与数值求解，只参与运动初始化，因此划分较大尺寸的网格，提高计算效率；
其次，确定鳍面的初始空间位置，并建立用户自定义函数(udf)链接到fluentmeshing中，用于控制鳍面初始化运动；
根据鳍面运动方程确定鳍面的初始空间位置，鳍面运动方程如下所示：
其中，参数d和s是鳍面上一点位置的空间坐标，θm是波动鳍运动幅值，t是波动鳍运动周期，φ是波动鳍运动相位差，λ是行波波长，无量纲系数ry和rz用于调节运动坐标；取时间t＝0s，根据方程可以确定波动鳍初始空间位置；将方程用c语言编写为用户自定义函数(udf)链接到fluent中，从而控制鳍面初始化运动；
最后，将得到的粗糙网格导入软件fluentmeshing，通过动网格技术进行网格运动并更新网格，得到粗糙网格仿真计算物理模型；
步骤4)，对粗糙网格仿真计算物理模型进行精细网格划分得到精细网格仿真计算物理模型；
如图4所示，将步骤3)得到的粗糙网格仿真计算物理模型进行格式转换符合网格划分软件格式，即把得到的粗糙网格仿真计算物理模型生成后缀名为cas的文件，导入软件fluentmeshing中进行精细网格划分；为了保证计算结果的准确性，精细网格总数不少于200万；将粗糙网格仿真计算物理模型划分为两个计算域，靠近仿生波动鳍模型的网格进行加密，网格类型为四面体网格，设置为一个计算域；远离仿生波动鳍模型的流场，网格类型为六面体网格，设置为另一个计算域；两个计算域通过共节点面连接；
步骤5)，对精细网格仿真计算物理模型进行仿真参数和运动条件设置；
将步骤4)得到的精细网格仿真计算物理模型导入软件fluentmeshing中，进行数值求解相关参数设置，采用非定常和压力基的湍流求解器求解流场；
设置的仿真参数包括涡粘系数模型、压力速度解耦算法、离散格式以及松弛因子；波动鳍运动条件设置指确定波状鳍运动方程，并编写相关用户自定义函数(udf)链接到fluentmeshing中，用于控制鳍面的周期运动；相关参数设置完毕后进行数值计算求解；
采用雷诺平均纳维-斯托克斯rans方程求解流场：
其中，
涡粘系数模型采用realizablek-εmodel，壁面函数采用yplus不敏感处理；压力速度解耦算法采用piso算法，时间项采用一阶隐式欧拉格式，压力梯度插值采用二阶格式，对流项采用二阶迎风格式；压力项的松弛因子降低到0.15以提高数值计算的收敛性；其中涡粘系数模型的运输方程为：
其中，gk和gb分别是由于流场平均速度和浮力导致的湍动能生成项，ym代表可压缩湍流中脉动扩张对总耗散率的贡献，σk和σε是湍流普朗特数，c2,c1ε和c3ε是常数，sk和sε是用户自定义的源项；
步骤6)、基于步骤5)设置的仿真参数和运动条件参数计算仿生波动鳍推力；
具体的，选定总计算时间、时间步长和计算迭代次数，通过软件fluentmeshing开展数值计算；计算过程中依据波动鳍推力的数学表达式，对推力变化曲线进行提取和监视；
具体的，如图5所示，基于步骤5)设置的仿真参数和运动条件参数，选定总计算时间为10s，计算步长0.005s，内迭代次数40次，采用迭代法进行仿生波动鳍推力计算，在数值计算求解过程中，采用分离式求解算法，在初值条件下先对流场速度进行求解，并进行压力矫正，使得求解的物理量满足守恒方程；其次对湍动能标量进行数值计算，在流场物理量求解完毕后，根据求解值对推力进行计算；波动鳍产生的推力由压力和摩擦力两部分组成；压力大小取决于仿生推进器的形状，而摩擦力则受流体粘度的影响，总推力可为：
ft＝a·fp a·ff
其中，fp是仿生推进器表面压力，ff是总摩擦力，a是方向向量；
由于波动鳍推力是一种效果力，虽然由上式定义，但不能直接求解，需要先求解压力和摩擦力；压力和摩擦力对波动鳍产生推进作用或者阻碍作用，总推力数值实际上由推进作用和阻碍作用两部分构成，即
ft＝t(t)-d(t)
其中t(t)和d(t)分别代表压力和摩擦力产生的推进作用和阻碍作用；
t(t)和d(t)分别为：
其中，p是单位压力，τj1是x方向的粘性应力张量分量，nj是方向向量，a是波动鳍总面积；根据上式，可以在求解过程中分别数值计算出t(t)和d(t)，从而计算出推力，并且监视推力曲线；
步骤7)，对仿生波动鳍推进性能进行分析；
推进性能指仿生推进器运动产生的推力；通过改变波动鳍运动的幅值与频率，获取波动鳍推力的变化规律；在一个周期t内，选取t＝0，t＝0.2t，t＝0.4t，t＝0.6t，t＝0.8t，t＝1t六个时刻的压力云图进行分析；选取t＝4t，t＝8t，t＝12t，t＝16t，t＝18t，t＝20t六个时刻的漩涡云图进行分析；记录不同频率下波动鳍产生的推力，绘制相应推力曲线。
绘制的波动鳍压力云图如图6所示，鳍面压力呈现周期变化规律；一个周期内鳍面高压和低压区域交替从鳍面脱落，高压区域出现在鳍面行波波谷处，低压区域出现在鳍面行波波峰处；一个周期内，一个完整的行波将从波动鳍鳍面向后脱落；由于波动鳍长度和行波波长在本发明中相同，因此鳍面始终存在一个完整的行波；当一个行波从鳍面脱落，会给周围流体施加一个向后作用力，从而流体给波动鳍施加一个向前作用力，产生向前的推力；
步骤7)中绘制的波动鳍周围漩涡云图如图7所示，波动鳍侧后方有明显漩涡从鳍面脱落，向后形成射流区域；射流区域的流体向后运动，从而推进仿生推进器向前运动；漩涡采用lambda-2法则进行识别；
步骤7)中绘制的波动鳍不同频率下平均推力曲线如图8所示，平均推力与波动频率近似成二次关系，根据流体力学中作动盘理论，该关系满足相应数学表达式。
综上所述，本发明提出的一种仿生波动鳍的推进性能数值计算方法，相比与传统数值计算方法，采用了两次网格划分方法，精确控制波动鳍初始位置和周期运动，高效计算了波动鳍推进性能，实现了对仿生推进器水动力性能的数值预测。
以上内容仅为说明本申请发明的技术思想，不能作为限定本发明的保护范围的依据，凡是按照本发明提出的设计构思和技术特征，在技术方案上所做的任何修改与替换，均在本发明权利要求书的保护范围之内。
1.一种仿生波动鳍的推进性能数值计算方法，其特征在于，包括以下步骤：
步骤1)，建立待计算仿生波动鳍的仿真计算物理模型；
步骤2)，确定仿真计算物理模型外流场范围以及边界条件；
步骤3)，根据外流场范围以及边界条件建立粗糙网格并进行波动鳍运动初始化，得到粗糙网格仿真计算物理模型；
步骤4)、对粗糙网格仿真计算物理模型进行精细网格划分得到精细网格仿真计算物理模型；
步骤5)，对精细网格仿真计算物理模型进行仿真参数和运动条件设置；
步骤6)、基于步骤5)设置的仿真参数和运动条件参数计算仿生波动鳍推力：选定总计算时间、时间步长和计算迭代次数，采用迭代法进行仿生波动鳍推力计算，得到最终仿生波动鳍推力性能数值。
2.根据权利要求1所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，确定外流场区域范围指确定波动鳍周围流场范围，采用尺寸大小确定波动鳍周围流场范围；波动鳍流场边界条件指流场前端半圆面设置为速度入口条件，后端平面设置为压力出口条件，周围圆柱面设置为壁面条件。
3.根据权利要求1所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，对仿真计算物理模型进行网格划分，确定鳍面的初始空间位置。
4.根据权利要求3所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，根据鳍面运动方程确定鳍面的初始空间位置，鳍面运动方程如下所示：
其中，参数d和s是鳍面上一点位置的空间坐标，θm是波动鳍运动幅值，t是波动鳍运动周期，φ是波动鳍运动相位差，λ是行波波长，无量纲系数ry和rz用于调节运动坐标；取时间t＝0s，根据方程可以确定波动鳍初始空间位置。
5.根据权利要求1所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，对步骤3)得到的粗糙网格仿真计算物理模型进行精细网格划分：将粗糙网格仿真计算物理模型划分为两个计算域，靠近仿生波动鳍模型的网格进行加密，网格类型为四面体网格，设置为一个计算域；远离仿生波动鳍模型的流场，网格类型为六面体网格，设置为另一个计算域；两个计算域通过共节点面连接，精细网格总数不少于200万。
6.根据权利要求1所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，将步骤4)得到的精细网格仿真计算物理模型导入软件fluentmeshing中，进行数值求解相关参数设置，采用非定常和压力基的湍流求解器求解流场。
7.根据权利要求6所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，采用雷诺平均纳维-斯托克斯rans方程求解流场：
其中，
8.根据权利要求6所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，涡粘系数模型采用realizablek-εmodel，壁面函数采用yplus不敏感处理；压力速度解耦算法采用piso算法，时间项采用一阶隐式欧拉格式，压力梯度插值采用二阶格式，对流项采用二阶迎风格式；压力项的松弛因子降低到0.15以提高数值计算的收敛性；其中涡粘系数模型的运输方程为：
其中，gk和gb分别是由于流场平均速度和浮力导致的湍动能生成项，ym代表可压缩湍流中脉动扩张对总耗散率的贡献，σk和σε是湍流普朗特数，c2,c1ε和c3ε是常数，sk和sε是用户自定义的源项。
9.根据权利要求1所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，采用分离式求解算法，在初值条件下先对流场速度进行求解，并进行压力矫正，使得求解的物理量满足守恒方程；其次对湍动能标量进行数值计算，在流场物理量求解完毕后，根据求解值对推力进行计算：ft＝t(t)-d(t)
其中t(t)和d(t)分别代表压力和摩擦力产生的推进作用和阻碍作用；
t(t)和d(t)分别为：
其中，p是单位压力，τj1是x方向的粘性应力张量分量，nj是方向向量，a是波动鳍总面积。
10.根据权利要求1所述的一种仿生波动鳍的推进性能数值计算方法，其特征在于，根据步骤6)得到的仿生推进器运动产生的推力对仿生波动鳍推进性能进行分析：通过改变波动鳍运动的幅值与频率，获取波动鳍推力的变化规律。
技术总结