基于CFD 水动力参数的水下机器人轨迹跟踪控制

罗一汉，吴家鸣，周汇锋

华南理工大学 土木与交通学院，广东 广州 510641

0 引 言

随着海洋资源的开发和水下领域作业任务增加，水下机器人扮演着越来越重要的角色。但机器人在水下环境工作时，有着非线性强、易受环境影响的特点。为获取良好的作业效果，对水下机器人的运动进行有效控制，使之能较好地跟踪预先设计的轨迹变得十分重要，这也使得轨迹跟踪控制成为水下机器人领域的主要研究热点之一。同时在对控制器控制效果进行检验时，水下机器人水动力参数的精确与否直接影响控制器控制效果的可信度。

针对水下机器人的水动力参数获取和航行轨迹跟踪控制，国内外有着相当丰富的研究成果。在水动力参数获取方面，使用CFD 技术获取是主要手段。Ramírez-macías 等[1]使用CFD 技术对某遥控水下机器人（ROV）艏部面罩进行了水动力模型建模，采用黏性流解算器ReFRESCO 求解单相湍流不可压缩流的稳态和非稳态N-S 方程，对不同来流方向下的ROV 横摇及其绕X，Y，Z轴旋转时的力和力矩进行了计算。程啸鹏等[2]基于Fluent 软件探究某小型水下机器人水动力与速度、角加速度之间的关系，并通过实体试验证明了仿真结果的可信性。赵凯凯等[3]以升力、阻力及纵倾力矩为对比指标进行数值仿真，对蛇形机器人外形进行优化，提高了此水下机器人的水动力性能。吴家鸣等[4]基于动网格和滑移网格技术，采用Fluent 软件对潜器做回转运动时在空泡影响下的螺旋桨推力特性进行了研究。陈东军等[5]基于CFD技术对带缆水下机器人系统升沉运动时的螺旋桨推力特性、水动力响应进行了分析研究。

在水下机器人轨迹跟踪方面，薛乃耀等[6]基于神经网络、滑模控制对轨迹跟踪提出一种变增益抗饱和系统并利用Matlab 进行仿真实验，结果表明，轨迹跟踪时推进系统推力饱和持续时间降低27%，从而验证了控制策略的可行性。Zhang等[7]提出一种新型自适应非线性二阶滑模控制器用于消除抖振运动，并使用自适应律实时估计外界干扰的幅值，实模试验结果表明，设计的控制器相较于线性二阶滑模面控制器有着更快的收敛速度和抗抖振能力。Liang 等[8]针对参数不确定性和外部干扰下的三维轨迹跟踪控制问题，提出了一种基于反步法和滑模的自适应鲁棒控制策略，使用模糊逻辑逼近非线性函数。仿真结果表明，该控制器有着良好的鲁棒性。陈伟等[9]针对欠驱动自主水下航行器基于视线法（LOS）和反步法设计了一种水平面轨迹跟踪控制器，并通过数值仿真验证了该控制器的有效性。

上述关于控制策略研究文献中的水动力参数来源有：CFD 计算获取[6]、实模参数[7-8]以及为验证控制器的有效性而适当选取[9]。然而，这些水动力参数存在精度差异，且设计控制参数数量过多也会导致调参困难的问题。鉴此，本文将使用CFD 软件STAR-CCM+和ANSYS AQWA 来获取所需的水动力参数，以保障参数的精度。在设计控制策略时，首先采取类似于文献[10-11]中使用的坐标变换方式，获取水下机器人水平面运动轨迹跟踪控制偏差系统。然后，为弥补文献[10]中设计的控制策略无法跟踪转艏运动期望角速度为0 时的运动轨迹缺陷，本文将基于反步法和滑模控制技术提出一种新的水下机器人水平面轨迹跟踪控制策略，并使用线性和非线性目标轨迹对控制器的有效性进行验证。

1 模型及参数获取

1.1 坐标系及运动方程

本文选用的坐标系如图1 所示。在固定坐标系O−XYZ中，规定OZ的正向指向地心，OX和OY轴在水平面内相互垂直，可选择任意方向为轴的正方向，选取水下机器人的主航向为OX的正向。运动坐标系o−xyz建立在机器人本体上，随机器人位置的变化而变化。为便于讨论，将坐标系原点定于载体重心处。ox轴与主对称轴一致，以艏部方向为正向；oy轴平行于基线面并垂直于ox轴，以指向右舷为正；oz轴与ox,oy轴垂直，指向底部为正。两个坐标系满足右手螺旋法则。

图1 固定坐标系与运动坐标系示意图Fig. 1 Sketch of fixed coordinate system and motion coordinate system

由给定的固定坐标系和运动坐标系，定义2 个变量 η和 υ，分别描述两个坐标系下机器人各自由度上的运动坐标值。两个坐标系的定义如下：

式中：x,y,z为机器人的位移；φ,θ,ψ为机器人的横摇角、纵摇角和艏摇角；u,v,w,p,q,r分别为各自由度上的运动线速度与角速度。二者之间的关系为

式中，R为两个坐标系的转换矩阵，定义如下：

其中，

由上式，水下机器人六自由度运动方程为

本文主要研究的是机器人水平面运动轨迹跟踪控制问题，仅考虑机器人在纵荡 (x)、横荡 (y)和转艏(Rz)这3 个自由度上的运动方程，其他自由度上的位移、速度及速度的一阶导数均取为0。选用的仿真模型近似三平面对称，在运动分析中部分参数矩阵仅考虑对角线部分，运动坐标系的原点为机器人的质心。简化后的机器人水平面运动方程为

式中，ηs={x,y,ψ}T；υs={u,v,r}T；阻尼力矩阵Ds=diag{D1,D2,D6}，其中Di=Dnoni·|υi|+Dlini，i=1,2,6；其他项也作相应简化。可见，若要对机器人运行轨迹进行控制，需要获取机器人的附加质量、非线性阻尼力项和线性阻尼力项的数据。

1.2 几何建模及其仿真数学模型水动力系数的确定

图2 为本文水下机器人的几何模型，机器人主体为一个艏艉近似半球形的过渡圆柱形外型结构（总长0.4 m，直径0.2 m），主体左右舷对称设置有一对推进导管螺旋桨，艏艉反对称设置有一对用于轨迹与姿态操纵的导管螺旋桨，分别用于提供进退运动的推力和产生转艏运动的力矩。螺旋桨动力由脐带缆提供，但假定脐带缆缆径较小，受流体干扰影响小，在运动过程中仅提供螺旋桨所需动力，对机器人运动的影响可忽略。研究中所需的机器人质心、惯性矩等相关参数根据图2的水下机器人几何模型从SolidWorks 建模软件中获取，具体数值如表1 所示。

图2 水下机器人模型Fig. 2 Underwater vehicle model

表1 水下机器人的相关参数Table 1 Relevant parameters of the underwater vehicle

为获取模型的水动力参数，本文选用CFD 软件STAR-CCM+和ANSYS AQWA 分别计算机器人模型阻尼力项参数和附加质量参数。在阻尼力项计算中，纵（横）荡的运动参数采用DFBI 模型计算，转艏运动阻尼力矩计算模型选用多参考系模型（MRF）计算。图3 为计算这些运动参数所使用的网格示意图。对于本文的纵荡运动阻尼力参数计算，计算域大小（x×y×z）取为4 m×3 m×2 m。为能够捕获流动细节，在水下机器人外部增设了2 层加密网格，网格类型为TRIM 类型。水下机器人主体网格类型为多面体网格，外部共有4 套网格，整体网格总数约为120 万（图3(a)）。计算域入口边界设置流速为0～1.5 m/s 的速度入口，出口边界设为压力出口，其他边界设为对称面，水下机器人壁面设为无滑移边界。计算模型选用Reliablek-ε湍流模型 ( RKE )。本文的横荡阻尼力计算方法与前述纵荡运动的一致。

转艏运动阻尼力参数计算所采用的计算域为高6 m、半径3 m 的圆柱体，轴线与z轴平行（图3(b)）。水下机器人外部共有4 层网格，与阻尼力项算例设计相似，另增设2 层TRIM 类型网格用于加密。机器人主体设定为旋转速度为0～2.5 rad/s的旋转域，网格类型为多面体网格，计算域边界设置参考阻尼力项算例，计算网格总数约120 万。

在计算过程中，将螺旋桨作为无转动浮体来考虑，故在计算机器人附加质量参数时采取了与文献[12]相同的方法对模型进行近似简化处理，以保证AQWA 计算效率和质量。图4 为所采用的模型网格示意图。附加质量参数计算采用ANSYS AQWA 软件，水密度1 025 kg/m3，计算域（x×y×z）为10 m×10 m×40 m，网格总数17 362。

图4 附加质量计算网格图Fig. 4 Grid meshes of additional mass calculation

1.3 网格无关性验证计算

为降低网格数对计算结果精度的影响，设计了少、中、多3 种网格数情况。对于阻尼力（矩）计算，网格总数分别取为 60 万、120 万和300 万，而附加质量计算的网格总数则分别取为0.8 万、1.7 万和2.6 万。通过对不同网格总数下的机器人运动阻尼力（矩）和附加质量计算结果的对比（图5），得到不同网格数下各水动力参数计算结果误差对比结果（表2）。表2 中纵（横）荡阻尼力和转艏阻尼力矩所得数据分别为流速0.9 m/s 和旋转速度1.5 rad/s 时的计算结果，附加质量项选取为附加质量项m1的计算结果，取网格数（多）情况下的计算结果作为对比标准。

由图5(a)和图5(b)可见，在阻尼力（矩）计算中，3 条曲线非常接近，结合表2 中阻尼力（矩）的误差对比，以纵荡阻尼力和转艏阻尼力矩为例，可见60 万和120 万网格数计算结果相较于300万网格数时的结果，纵荡阻尼力项相对误差由2.314%减少至0.833%，转艏阻尼力矩项由2.271%减少至0.189%，这表明在同一设计速度点下，3 种网格数下获取的阻尼力（矩）值相差不大。因此，虽然增加网格数可获得更精确的计算结果，但收益不大。由图5(c)和表2 还可见，对于附加质量计算结果，3 个参数在3 种网格数下的值基本保持不变，相对误差在2%左右。因此，综合考虑计算资源和计算精度，本文阻尼力项和附加质量项选择为120 万网格数和1.7 万网格数的计算结果。此外，本文对阻尼力项中的非线性项仅考虑了二次非线性项，对所得阻尼力（矩）参数进行二次多项式拟合后所得到如图6 所示拟合曲线。进而可得阻尼力矩阵为：Dnons=diag{12.331·|u|,21.417·|v|,0.612·|r|}，Dlins=diag{1.012,1.541,0.137}，附加质量矩阵为：MAs=diag{3.788,15.302,0.382}。拟合曲线表明二阶非线性阻尼项远大于线性阻尼项，只有在机器人运动速度很小的情况下，才能忽略非线性项阻尼对机器人运动的影响。

图5 不同网格数下所得水动力参数对比Fig. 5 Comparison of hydrodynamic parameters under different grid numbers

表2 不同网格数下水动力参数计算结果对比表Table 2 Comparison of calculated hydrodynamic parameters under different grid numbers

图6 阻尼力（矩）拟合曲线示意图Fig. 6 Diagram of damping force (moment) fitting curve

2 水平面轨迹控制策略

通过CFD 技术获取仿真模型的水动力参数后，为使用模型对机器人水平面轨迹跟踪进行控制，还需对作用力矩阵Fs进行定义。本文研究的水下机器人配置有4 个螺旋桨，其中主体左右舷的两螺旋桨仅提供x方向上的推力，艏艉端的2 个螺旋桨反对称布置，为转艏运动(Rz)提供一个力偶矩。因忽略脐带缆对机器人运动的影响，且y自由度上无作用力，故本系统为一欠驱动系统。此时，作用力矩阵可表示为

式中：T为x方向上的推力；N为转艏力矩。

结合式 (7)和式 (8)，可得如下运动系统：

式中，符号定义同式 (7)，并令x={x,y,ψ,u,v,r}T。参考文献[10-11]采用的坐标变换形式，可得：

式中，z6d为转艏运动角速度的期望值，在工程实际中该值一般变化速度不太大，可近似认为z˙6d≈0。因此，第1 个式子对时间求导，将剩余两式代入，可得：

3 水下机器人轨迹跟踪控制特性分析

本节在Matlab/Simulink 平台对两类轨迹（前者为线性轨迹，后者为非线性轨迹）下控制策略的有效性进行验证。具体参数分别为：

3.1 线性轨迹跟踪

已知式(27)所示线性轨迹仿真目标参数及式(29)和式(30)所示水下机器人系统参数，且控制参数选择为：c1=10,η1=5,c2=15,η2=4,k=5,γ=10。但初始条件中z6d=rd=0，若依据前述控制策略分析会存在一定的跟踪偏差，故保持其他参数不变，选取z6d=rd=0.6sin(0.1t)作为对照组进行仿真计算，结果如图7～图8 所示。

由图7 可见，在设置的不同z6d条件下，e1，e3，e6响应总是能较快收敛到0，时间约为2.2 s；e4和e5在两种情况下也都能收敛到0，而e2在z6d=0的条件下不为0，但会收敛至一个稳定值，而当z6d≠0时，可收敛至0，这与理论分析相符。由图8(a)也可见，在z6d=0的条件下，仿真轨迹与目标轨迹存在一定的偏差，但两类轨迹相互平行，表明虽然在x，y自由度上存在一定的位置偏差，但艏摇角ψ可被跟踪上，即e3=0；而当z6d≠0时，仿真轨迹在一段时间后会与目标轨迹重合，实现了较好的追踪效果。由图8(b) 所示的螺旋桨推力示意图可明显看出：2 种条件对左右舷提供进退运动推力的螺旋桨影响不大，数值稳定在7 N 左右，但对机器人艏艉端提供偏航力矩的两个螺旋桨推力的作用较大，可见螺旋桨所需提供的推力明显提升，而变化较为平稳。因此，本文所提控制策略能实现对线性运动轨迹的跟踪，可解决文献[10]中轨迹跟踪控制器无法跟踪线性轨迹（z6d=0）的问题。

图7 不同z6d 条件下偏差项响应图Fig. 7 Response of deviation under different z6d conditions

图8 不同z6d 条件下运动轨迹和螺旋桨推力响应图Fig. 8 Response of trajectory and propeller thrust under different z6d conditions

3.2 非线性轨迹跟踪

非线性轨迹仿真目标参数如式(28)所示，水下机器人系统参数如式(29)和式(30)所示，控制参数选择为：c1=2,η1=2,c2=1.5,η2=2.5,k=5,γ=10，仿真结果如图9～图12 所示。

从图9 可得，在非线性轨迹跟踪过程中，误差项ei(i=1 ∼6)，都能收敛到0，在10 s 左右便收敛至0 附近，表明控制策略的控制效果较好。但同时可见部分偏差参数存在较大超调，结合图10 轨迹跟踪图可见，由于在初始时刻仿真轨迹和目标轨迹存在较大偏差，使得在偏差收敛过程中有较大的起伏变动，而当偏差收敛至0 后，曲线变得平稳。进一步观察图10，可见仿真轨迹能够快速跟踪上目标轨迹，并保持了良好的持续跟踪效果，趋近过程中的最大超调量约10%。由图11 可更直观地看出，在初始时刻实际运动轨迹和目标轨迹的偏差最大，这与所设计的目标轨迹和机器人初始位置状态相符。在运动过程中，伴随着控制策略的控制动作，约在10 s 左右使得两个轨迹位置偏差达到0，且后续几乎无变化，表明此时已较好地实现了轨迹跟踪，并能保持持续的跟踪效果。由图12 所示的螺旋桨推力示意图可见：在0～10 s 内，左右舷螺旋桨推力存在明显的调节过程，而在偏差项都收敛至0 后，螺旋桨推力变化较为平缓，艏艉端螺旋桨的推力在约3 s内便趋于稳定，这与图9 所示的偏差项变化趋势是相符的。系统稳定后，左右舷螺旋桨推力维持在3.2 N 左右，艏艉端螺旋桨推力维持在15 N左右。因此，对于非线性轨迹，本文设计的控制策略依旧能够实现较好的追踪效果并保持螺旋桨推力稳定。

图9 偏差项响应图Fig. 9 Response of deviation term

图10 轨迹跟踪示意图Fig. 10 Diagram of trajectory tracking

图11 轨迹跟踪位置偏差示意图Fig. 11 Diagram of tracking trajectory deviation

图12 螺旋桨推力示意图Fig. 12 Diagram of propeller thrust

综合上述两类运动轨迹的仿真结果，表明本文设计的控制器可使水下机器人快速跟踪目标轨迹，能够保持较好的持续跟踪效果，验证了控制器的有效性；两种轨迹下的螺旋桨推力在系统稳定后都趋于稳定，但线性轨迹下的螺旋桨推力存在峰值，而非线性轨迹下的推力较为平缓，这可能是因为初始时刻所设计的线性轨迹变化速度较快，使得螺旋桨必须提供一个较大推力以满足轨迹快速趋近的要求所致。因此，在设计实际轨迹时，应充分考虑推进器作业能力的大小。

4 结 论

本文针对水下机器人水平面轨迹跟踪问题进行了研究。首先，考虑二阶非线性阻尼力推导了水下机器人三自由度水动力数学模型，使用STAR-CCM+和ANSYS AQWA 软件分别获取了所需阻尼力和附加质量参数。然后，结合模型自身的动力配置条件、计算所得的水动力参数，推导得到了机器人水平面运动轨迹跟踪控制偏差系统。其后，将偏差系统分为3 个子系统，基于反步法和滑模控制技术逐个进行控制策略设计，提出了一种新的控制策略，使得能够实现对线性运动轨迹（转艏运动期望角速度为0）和非线性运动轨迹的跟踪。最后，使用设计的控制器，针对上述两类运动轨迹，在Matlab/Simulink 平台中进行了轨迹跟踪仿真计算。得到主要结论如下：

1) 在水动力参数计算方面，本文所得阻尼项拟合曲线表明，二阶非线性阻尼项远大于线性阻尼项，只有在机器人运动速度很小的情况下才能忽略二阶非线性阻尼项对运动的影响。

2) 对于线性运动轨迹（转艏运动期望角速度为0）跟踪问题，本文设计的控制策略可弥补所对比控制策略无法跟踪线性运动轨迹的缺陷。对于非线性轨迹跟踪问题，本文设计的控制策略依然可实现较好的跟踪效果，趋近过程中的最大超调量约为10%，轨迹跟踪所需时间约为10 s。

3) 在两类运动轨迹下，螺旋桨推力在系统稳定后都能趋于一个稳定值，但线性轨迹时的螺旋桨推力曲线存在峰值，而非线性轨迹时螺旋桨推力曲线较为平缓，其原因是本文设计的轨迹参数变化较快所致。因此，在设计实际运动轨迹时，应充分考虑推进器作业能力的大小。