深海运载器无动力纵倾上浮运动特性研究

赵志超,李天辰,谷海涛,高浩

(1.中国科学院 沈阳自动化研究所机器人学国家重点实验室,辽宁 沈阳,110016;2.中国科学院 机器人与智能制造创新研究院,辽宁 沈阳,110169;3.中国科学院大学,北京,100049;4.海军研究院,北京,100161)

0 引言

深海运载器是一种将无人机等任务载荷从深海环境快速投送到指定海域的运载工具。运载器携带任务载荷时,采用无动力上浮方式可有效节省能源消耗,但由于缺少动力装置推进,仅依靠压载和浮力调节产生的有限净浮力驱动上浮,存在大深度上浮耗时长的特点。无动力上浮按运动形式可分为无纵倾上浮、纵倾上浮和空间螺旋上浮。通常水下运载器采用无纵倾上浮和空间螺旋上浮比纵倾上浮产生的垂向速度小,如蛟龙号载人潜水器通过调整上浮压载水舱,可实现无纵倾上浮,最大上浮垂向速度约0.813 m/s[1];水下滑翔机最大滑翔速度不超过1.3 m/s,而空间螺旋上浮垂向速度不超过0.5 m/s[2-3]。无纵倾上浮、小纵倾上浮或空间螺旋上浮形式无法保证深海无动力运载器快速上浮8 kn 以上的速度要求。为此,文中研究以大纵倾角或垂直姿态上浮方式,实现运载器快速上浮垂向速度的要求。

深海运载器大深度上浮过程中的水动力和环境扰动复杂,增大了运载器无动力上浮运动预报难度。目前,水下运载器上浮运动预报研究主要采用模型试验法、计算流体力学(computational fluid dynamics,CFD)计算法以及数学模型方法[4]。其中模型试验法和CFD 数值计算方法较为成熟。张天健等[4]采用数值计算和快速预报方法,研究了大深度浮力驱动式运载器上浮运动规律。张海洋等[5]采用数值计算方法,研究了无动力运载器自身比重、重浮心位置等对纵倾上浮状态的影响。由于大深度上浮运动预报受模型试验条件和CFD 数值计算规模限制,无法快速完整地分析大深度上浮过程浮力变化和洋流扰动等对运载器实时和总体上浮状态的影响,因此相关研究逐渐采用数学模型方法。采用数学模型方法进行深海运载器的大深度纵倾上浮运动预报的主要难点涉及:欧拉角描述的空间运动方程在大纵倾角或垂直上浮时姿态求解奇异问题、大深度上浮过程中浮力变化和洋流干扰的数学建模[6]。针对上述问题已有相关研究:徐正武等[6]将四元数法应用于水下航行器的运动建模,有效避免了使用欧拉法描述空间运动的奇异问题。王海涛等[7]建立了随机海浪扰动模型,并使用四元数法解决了运载器垂直出水姿态解算奇异问题,研究了不同海况与主浪向角对运载器垂直出水状态的影响;耿斌斌等[8]基于大量水动力试验数据,构建了无动力水下运载器的六自由度数学仿真模型,讨论了不同初始状态下偏转、速度变化等因素对运载器运动状态的影响;高伟等[9]基于自主水下航行器(autonomous undersea vehicle,AUV)空间运动方程,建立了稳态无动力下潜运动的仿真方法,揭示了AUV 各状态量与控制量及海洋环境之间的变化规律。但上述类似研究未探究无动力运载器采用大纵倾角或垂直姿态进行大深度上浮过程中,洋流、浮力变化和发射条件影响下的上浮运动特性。

为此,文中以某型回转体AUV 为计算模型,基于四元数法和AUV 空间运动方程,构建了运载器空间运动仿真模型,并添加大深度上浮过程洋流和浮力变化干扰,在该模型基础上探究了净浮力、重浮心距离、舵角、发射条件以及洋流和浮力变化扰动下,运载器做大纵倾角和垂直姿态上浮时的运动特性。

1 数学模型

1.1 坐标系定义

为描述深海运载器纵倾上浮运动,可引入惯性坐标系E-ξηζ和体坐标系G-xyz,惯性坐标系原点位于海平面,体坐标系原点位于运载器浮心位置,各坐标轴均按右手系确定,如图1 所示。运载器的六自由度运动可通过以下向量描述:惯性系中运载器浮心位置坐标η1=[ξ η ζ]T,速度VE=[VξVηVζ]T,姿态角η2=[φ θ ψ]T,体坐标系下浮心速度V=[u v w]T,角速度 Ω=[p q r]T,运载器所受合外力F=[X Y Z]T,合外力矩M=[K M N]T。

图1 深海运载器坐标系示意图Fig.1 Coordinate system of a deep-sea vehicle

1.2 四元数法

文中研究的深海运载器采用大纵倾角或垂直姿态上浮,即上浮纵倾角接近或超过90°。运用四元数法描述水下航行器的姿态,能够克服欧拉角描述的空间运动方程奇异性的缺点[10]。四元数是Hamilton 于1843 年提出的一种超复数。20 世纪70 年代四元数开始在控制工程中得到应用。四元数包含一个三维矢量和一个标量。根据欧拉旋转定理,刚体绕固定点的任意位移可以绕通过此点的某一轴转动某一角度而得到。四元数的矢量部分表示转动轴的方向,其标量部分则与绕转动轴转过的角度有关,因此定义单位四元数为

式中:λ为转动轴的单位矢量;β为绕转动轴转过的角度。

1.3 基于四元数法的运载器空间运动方程

使用四元数法代替欧拉角描述运载器浮心运动的关系式为

描述运载器浮心旋转的运动学方程,通过四元数微分方程(4)实现,其中E2(ε)如式(5)[10],即

将式(2)和式(4)表示为矩阵形式,得到基于四元数的运载器空间运动学方程

深海运载器在无动力上浮过程中,所受的合外力(矩)包括:水静力(矩)、水动力(矩)以及舵的操纵力(矩)。根据文献[10-12],深海运载器的动力学方程可表示为

式中:M为深海运载器质量、转动惯量、惯性水动力系数和重浮心距离组成的系数矩阵,如式(8);Fhs和Mhs为运载器重力和浮力作用产生的水静力和力矩[10];Fhd和Mhd为运载器受到的粘性水动力和力矩;Fop和Mop为运载器操纵舵产生的力和力矩;Fτ和Mτ为浮力变化等外部环境扰动产生的力和力矩。

式中:P和B分别为航行器所受到的重力和浮力;xg、yg和zg为运载器重心在体坐标系中的位置;xb、yb和zb为浮心在体坐标系中的位置。

粘性水动力和力矩可表示为[11]

1.4 四元数姿态解算

任意时刻运载器姿态都是确定的,因此欧拉角和四元数方法表示的姿态参数可以相互转换。运载器初始姿态定义和实时姿态解算过程中常用到欧拉角与四元数之间的转换关系。欧拉角表示的动坐标系向惯性坐标系的坐标变换矩阵[10]为

由于任意时刻四元数和欧拉角表示运载器姿态的坐标变化矩阵等价,对比式(3)和式(18)中的第1 列和第3 行元素可知

化简得四元数到欧拉角转换关系为

运载器的初始姿态可通过初始四元数定义,初始姿态角φ,θ,ψ到初始四元数σ1,σ2,σ3,η的转换关系为[13]

1.5 海洋环境扰动模型

1.5.1 洋流扰动模型

洋流的周期性约为几个月的量级,因此可以认为在给定的海况中,流速在大小和方向上不随时间变化,洋流的强度和方向只在不同深度的水层中会发生变化[14]。深海运载器大深度上浮过程中,若不考虑涡流和海流垂向的干扰,洋流速度随深度变化的规律[15]如图2 所示。洋流干扰可通过体坐标系中的相对速度Vr=[urvrwr]T施加到运载器本体上。假设惯性坐标系中的洋流速度Vc=[ucvcwc]T且=0,则运载器在体坐标系中相对洋流的速度[10]

图2 洋流速度随深度变化剖面图Fig.2 Ocean current velocity at different depths

1.5.2 浮力变化模型

深海运载器上浮过程浮力变化由载体排水体积和工作环境的海水密度决定[16]。运载器上浮过程周围海水温度、盐度和压力等参数随深度不断改变,引起运载器排水体积和海水密度变化,进而产生浮力改变。浮力变化计算依据遥控水下航行器(remote operated vehicle,ROV)Okeanos Explorer于2016 年6 月对马里亚纳海沟探测得到的6000 m温盐深传感器(conductivity,temperature,depth,CTD)数据[17],绘制海水温度和盐度随深度的变化关系,如图3 所示。海水的温度T和盐度S均在0～1000 m内变化剧烈,之后趋于平缓。盐度总体变化量较小,最大变化量约为1.13‰。

图3 海水温度和盐度随深度变化曲线Fig.3 Temperature and salinity values at different depths

假设深度ζ处海水的平均密度为ρζ,可用海水平均密度ρa近似代替。不考虑温度和盐度变化,则海水的压强变化量为

海水温度、盐度和压力对运载器上浮过程浮力变化ΔB的综合影响可表示为[18]

式中:ΔρT、ΔVT、ΔρP、ΔVp分别表示运载器上浮过程温度和压力变化引起的海水密度和排水体积的变化量;ΔS为盐度变化量。ROV 的耐压舱和结构件的材料几乎全是铝合金[16],这里将运载器作为同一材质考虑,浮力变化计算的相关参数,如表1[16,18]所示。

表1 运载器浮力变化计算参数Table 1 Calculation parameters of the buoyancy change of the vehicle

由表1 和式(16)可得,深海运载器上浮过程所受浮力变化如图4 所示,上浮过程深度从6000～644.2 m,浮力变化量为正,运载器所受净浮力增加了4.478 N;从深度644.2 m 上浮至水面过程中,浮力变化量为负,运载器所受净浮力减小了1.979 N,上浮过程运载器所受浮力增加约2.5 N。深度在1500～6000 m 区间内,运载器所受浮力变化量与深度呈正相关。

图4 浮力随深度变化曲线Fig.4 Buoyancy change at different depths

1.6 仿真建模

REMUS 100 AUV 主体是采用Myring 艇体轮廓方程的回转体外形,艉部十字形布置安装有4 个相同的NACA 0012 舵,其动力学模型准确性已通过多次海上试验验证,可作为文中研究的计算模型,其主要参数[11]:特征长度1332.7 mm,回转体最大直径191 mm,质量30.479 kg,全排水体积约0.0315 m3,总浮力约317 N,浮心距离艏部顶点0.611m,转动惯量为Ixx=0.177 kg·m2、Iyy=3.45 kg·m2、Izz=3.45 kg·m2。基于四元数法、AUV空间运动方程和海洋环境扰动模型,使用Simulink和S 函数建立了深海运载器空间运动仿真模型,模型结构及各模块组成如图5 所示。

图5 深海运载器空间运动仿真模型Fig.5 Space motion simulation model of the deep-sea vehicle

2 无动力纵倾上浮运动仿真分析

2.1 模型验证

为验证深海运载器空间运动仿真模型对大纵倾角上浮运动预报的准确性,基于STAR-CCM+及其动态流体相互作用(dynamic fluid body interaction,DFBI)数值计算方法,对运载器大纵倾角上浮运动进行动态仿真。CFD 模型设置平均海水密度ρ=1024 kg/m3,动力黏度µ=8.887×10-4Pa·s,湍流模型采用可实现的k-ε两层模型,入口速度和出口压力均设为0,运载器表面为无滑移壁面,并采用棱柱层网格进行划分,边界条件如图6(a)所示。运载器表面网格棱柱层数15 层,棱柱层延伸系数1.501,棱柱层厚度0.01 m,最大网格单元尺寸0.2 m,面网格最小尺寸0.05 m,面网格增长率1.3,最终得到网格单元总数92.283 万,面网格数量274.826 万,网格划分场景如图6(b)所示。DFBI设置中运载器六自由度体初始参数配置为重浮心纵向距离xgb=0.15 m,径向距离ygb=0.0196 m,净浮力W=18.069 N,深度ζ=6000 m,初速度u0=1 m/s垂直姿态发射,稳态垂向速度云图如图6(c)所示。

图6 基于STAR-CCM+的无动力纵倾上浮运动仿真Fig.6 Simulation of trim ascent motion of unpowered vehicle based on STAR-CCM+

相同工况下,使用基于Simulink 的空间运动仿真模型对运载器上浮运动进行预报,并与STARCCM+数值计算方法所得结果对比,如图7 所示,2 种方法计算所得稳态上浮垂向速度分别为-3.732 m/s 和-3.633 m/s,稳态纵倾角分别为81.085°和81.975°。上浮垂向速度和深度随时间的变化曲线基本吻合,稳态纵倾角计算结果相对误差小于1°,因此基于Simulink 的空间运动仿真模型,其计算精度满足工程应用要求。

图7 深海运载器纵倾上浮运动仿真对比验证Fig.7 Comparison verification of trim ascent motion of the deep-sea vehicle

为验证CFD 数值计算方法本身的准确性,配置运载器重浮心纵向距离xgb=0.15 m,径向距离ygb=0,净浮力W分别为0.05V,0.10V,0.15V,深度ζ=6000 m,初速度u0=1 m/s 垂直姿态发射。当运载器垂直上浮阻力等于净浮力时,上浮垂向速度达到最大,之后以最大速度保持稳态匀速运动,根据式(25)[11]可得不同净浮力条件下的模型阻力系数及其平均值,如表2 所示。已知REMUS 100 AUV拖曳水池试验测得的模型平均阻力系数Cd=0.12[11],而数值计算方法得到的阻力系数平均值与水池试验测得的结果误差较小,仅为1.873%,至此该数值计算方法本身的准确性得到了验证。

表2 不同净浮力下运载器阻力系数Table 2 The axial drag coefficient of vehicle at different net buoyancy

式中:W为净浮力;Af为最大投影面积。

综上所述,基于Simulink 建立的深海运载器空间运动仿真模型计算精度满足工程应用要求,具有较高的准确性,能够用于运载器大纵倾角和垂直上浮运动的特性研究。

2.2 控制量对纵倾上浮运动的影响

深海运载器纵倾上浮是指通过配置一定上浮参数,运载器上浮运动达到稳态后保持大纵倾角,以实现较大上浮垂向速度的运动形式。首先探究深海运载器在无海洋环境干扰下的纵倾上浮运动特性。采用控制变量法,通过调整运载器的重浮心纵向距离xgb、径向距离ygb、净浮力W、水平舵角δs、初始速度大小u0及初始纵倾角θ0,共6 个控制量,分析各控制量对运载器纵倾上浮状态参数的作用规律。

2.2.1 净浮力与重浮心纵向距离

净浮力和重浮心纵向距离是实现运载器纵倾上浮的主要控制量。忽略海洋环境扰动,纵倾上浮运动过程中保持运载器δs=10°和ygb=0 不变,以初始速度u0=1 m/s 垂直姿态发射,通过调整运载器的净浮力W=0～0.15V和重浮心纵向距离xgb=0～0.25L中的一系列固定值,计算运载器上浮运动达到稳态后的状态参数,探究运载器W和xgb对稳态上浮时纵倾角θ和垂向速度Vζ的作用规律,如图8 和图9 所示。

图8 稳态纵倾角与净浮力、重浮心纵向距离的关系Fig.8 Relationship among steady trim angle,net buoyancy,longitudinal distance between center of gravity and buoyancy

图9 稳态垂向速度与净浮力、重浮心纵向距离的关系Fig.9 Relationship among steady vertical velocity,net buoyancy,longitudinal distance between center of gravity and buoyancy

净浮力和重浮心纵向距离对稳态纵倾角的影响,主要通过水静力矩,粘性和惯性力矩及操纵力矩相互作用产生。由图8 可知,重浮心纵向距离一定时,稳态纵倾角随着净浮力的增大而减小。这是因为模型净浮力的调整是通过保持浮力不变减小重力实现的。假设运载器舵角不变,以某一稳态纵倾角稳态运动,此时通过减小重力增大净浮力的方式会导致水静力瞬时回复力矩减小,运载器在瞬时水动力作用下纵倾角逐渐减小,最终运载器在舵的操纵力矩、水静力回复力矩和水动力纵倾力矩平衡时以更小的稳态纵倾角实现了新的稳态运动。净浮力一定时,稳态纵倾角随着重浮心纵向距离的增大而增大。净浮力配置越小时,重浮心纵向距离对稳态纵倾角的作用范围越大,更易通过改变重浮心纵向距离,调整运载器的上浮姿态。净浮力越大时,运载器稳态纵倾角作用范围越小,此时重浮心纵向距离对稳态纵倾角调节范围越有限。

为实现运载器大纵倾角快速上浮,减小上浮过程垂向阻力,在有限净浮力条件下,可配置较大的重浮心纵向距离,以增大稳态上浮纵倾角。根据图9,当采用较大重浮心纵向距离时,净浮力与稳态上浮垂向速度存在二次关系。运载器固定水平舵角10°,当采用净浮力0.10V、重浮心纵向距离0.25L时,稳态上浮垂向速度可达到-3.598 m/s。综合考虑,对于大纵倾角快速上浮运动,可将运载器配置为净浮力0.10V,重浮心纵向距离0.25L,以实现较大的上浮垂向速度。

2.2.2 重浮心径向距离与水平舵角

忽略海洋环境干扰,运载器上浮过程保持δr、W=0.10V,xgb=0.25L不变,通过配置重浮心径向距离ygb分别为0.05D、0.15D和0.25D,水平舵角δs依次取-20°～+20°内固定值,计算运载器上浮运动达到稳态后的状态参数,探究重浮心径向距离和水平舵角对运载器稳态上浮参数的作用规律,如图10 所示。

由图10(a)可知,水平舵角一定时,重浮心径向距离对稳态纵倾角的影响较小,随着运载器重浮心径向距离增加,稳态纵倾角呈小幅减小的趋势。重浮心径向距离固定不变,稳态纵倾角随正向舵角的减小而增大,随负向舵角的减小而减小。若运载器配置重浮心纵向距离0.25L、径向距离0D,通过调整水平舵角为-20°～20°,可控制运载器稳态纵倾角在40°～140°之间。

重浮心径向距离一定时,稳态上浮垂向速度与水平舵角之间具有二次关系,如图10(b)。水平舵角从-20°减小到0°,后从0°增加到+20°过程中,稳态垂向上浮速度先增大后减小,并在0°舵角附近达到最大。最大稳态上浮垂向速度对应的水平舵角,随重浮心径向距离的增大而增大。

图10 稳态上浮参数与重浮心径向距离、水平舵角的关系Fig.10 Relationship among steady parameters,radial distance between center of gravity and buoyancy,rudder angle

当运载器配置净浮力0.1V、重浮心纵向距离0.25L、径向距离0D、水平舵角0°,运载器达到稳态上浮时的纵倾角为90°的垂直姿态,此时稳态上浮垂向速度达到最大-4.424 m/s。

2.3 环境扰动对纵倾上浮运动的影响

2.3.1 洋流扰动

忽略洋流对运载器的垂向干扰,将洋流扰动简化为均匀水平流,研究洋流扰动对运载器上浮运动轨迹的影响。运载器参数配置为W=0.1V、xgb=0.25L、ygb=0.1D,δs=0°,深度ζ=6000 m,以初速度u0=1 m/s 垂直姿态发射。对运载器本体施加不同方向的洋流速度扰动,得到运载器在不同方向来流下的浮心运动轨迹,Eζ轴数据取负值,如图11 所示。

图11 洋流影响下的运载器浮心运动轨迹Fig.11 Motion trajectory of vehicle under ocean current

对比不同来流干扰下,运载器在大纵倾角上浮运动过程中产生的水平漂移距离,如表3 所示。无偏顺流时,水平面纵向漂移距离最远;无偏逆流作用下,水平面纵向漂移距离最小。横向洋流扰动即顺流90°或逆流90°时,横向漂移距离与无洋流干扰时的纵向漂移距离相等;斜向45°顺流纵向距离相比斜向45°逆流时偏大,而横向漂移距离相等。

表3 洋流扰动下的深海运载器水平面漂移距离Table 3 Horizontal drift distance of deep-sea vehicles under the disturbance of ocean current

常见无动力运载器做大深度上浮运动的水平面漂移距离如表4 所示。其中载人潜水器采用空间螺旋下潜方式的水平漂移距离最小,但垂向速度受到极大限制;相对于采用小纵倾角上浮方式的水下滑翔机和深海AUV,深海运载器采用大纵倾角或垂直姿态上浮方式在保证快速上浮垂向速度要求的同时,能有效减少其水平面漂移距离。

表4 水下航行器纵倾上浮水平漂移距离Table 4 Horizontal drift distance of trim ascent of the undersea vehicle

2.3.2 浮力变化扰动

深海运载器大深度上浮过程中水文参数不断变化,所受浮力随深度变化如图4 所示。为研究浮力变化扰动特性,忽略洋流和初始发射条件对上浮状态参数的影响。保持运载器初始参数配置W=0.1V、xgb=0.25L、ygb=0.1D,δs=0°不变,深度ζ=6000 m,以初速度1 m/s 垂直姿态发射。对比分析有无浮力变化对运载器上浮过程纵倾角、垂向速度、水平速度等状态量的影响,如图12 所示。

图12 浮力变化对运载器纵倾上浮运动的影响Fig.12 The effect of buoyancy change on the trim ascent motion of vehicle

图12 中负号表示无浮力变化,正号表示考虑浮力变化影响,结合图4 可知,深海运载器能够依据浮力变化对上浮运动状态不断调整,上浮过程运载器纵倾角随浮力变化量的减小而增大,垂向和水平速度随浮力变化量的减小而减小。上浮过程前20 min,运载器所处深度在626～6000 m,浮力变化量为正,纵倾角相比无浮力变化时偏小,水平和垂向速度相比无浮力变化时偏大。20 min 后,浮力变化量为负,纵倾角相比无浮力变化时偏大,水平和垂向速度相比无浮力变化时偏小。但由于上浮过程浮力整体变化量为正,运载器所受净浮力总体增大,平均速度增加,因此上浮过程总耗时比无浮力变化时偏少。

2.3.3 初始发射条件

深海运载器的发射过程位于6000 m 深度海底,海底地形因素复杂多变。考虑海底地形导致运载器初始俯仰角偏差在±10°内,即初始俯仰角θ0=80°～100°,初始发射速度u0依次取1,3 和5 m/s,纵倾上浮过程保持W=0.1V、xgb=0.25L、ygb=0.1D、δs=δr=0°不变。忽略洋流和浮力变化干扰,探究初始俯仰角和发射速度对运载器纵倾上浮状态的作用规律,如图13 所示。

根据图13(a)和(b),初始速度大小和初始俯仰角不影响运载器稳态运动后的上浮状态参数。同一初始俯仰角发射,初始速度越大,纵倾角达到稳定前的摆动幅值越大,纵倾角稳定所需的时间就越长,但不影响运载器稳态上浮时的垂向速度和纵倾角。同一初始速度大小发射,初始俯仰角不影响运载器的稳态上浮垂向速度和纵倾角。

图13 初始发射条件对纵倾上浮运动的影响Fig.13 The effect of initial lauch conditions on the trim ascent motion

由图13(b)和(d)可知,运载器以初始俯仰角θ0=80°～100°发射,其上浮过程存在纵倾角等于或超过90°的情况,仿真模型仍能够准确地对运载器的上浮状态进行预报,验证了基于四元数法的深海运载器空间运动仿真模型,在运载器采用大纵倾角或垂直姿态上浮的运动预报过程中,能有效避免欧拉角表示的空间运动方程姿态求解奇异问题,提高了运载器在上浮过程中纵倾角等于或超过90°时的运动预报精度。

3 结束语

基于四元数法、AUV 空间运动方程及大深度洋流和浮力变化扰动等数学模型,通过MATLAB/Simulink 建立了深海运载器空间运动仿真模型,并通过CFD 数值计算方法验证了该模型的准确性。基于该模型设计了一系列仿真工况,探究了上浮控制参数和海洋环境扰动对深海无动力运载器大纵倾角和垂直上浮状态的作用规律。通过探究不同初始发射条件下运载器上浮状态,验证了四元数法在解决运载器大纵倾角和垂直上浮姿态求解奇异问题的有效性。相关研究为最终确定深海运载器的上浮形式、运动预报方法和总体布局提供了参考。后续需要进一步考虑洋流垂向扰动、波浪等环境因素对深海无动力运载器大深度纵倾和垂直上浮运动的影响。