0 引　言

台风引发的经济损失是全部自然灾害经济损失的重要组成部分，其引发的洪水灾害对沿海地区实际GDP、居民收入、消费等宏观经济指标均有负面影响^[1]。台风占我国自然灾害经济损失的1/3，由于预报能力不足，年均经济损失约1000亿元^{[2 - 3]}。目前关于台风预报的研究主要集中在大气科学领域，对于台风与海洋之间多尺度响应机制和反馈作用缺乏认识，导致台风路径和强度预报水平提升有限。在路径预报方面，Landsea等^[4]和Zhou等^[5]的研究表明，台风路径预报误差的减小趋势在近几年明显减缓。Yu等^[6]最新研究成果表明，假设过去10年的预报技术发展趋势能继续保持，预计将在2035年前后达到当前台风路径预报极限。在强度预报方面，缺乏台风过境时的现场有效观测资料，导致台风强度预报水平近几十年几无提升^[3]。2021年，《Science》发布“全世界最前沿的125个科学问题”，其中“人类能否更准确地预测灾害性事件（海啸、台风、地震）”被列为125个最前沿科学问题之一^[7]。科学家已经证实台风核心区海气界面海表飞沫与混合层水体热容量是影响台风路径和强度预报的重要因素^{[8 - 9]}。

传统海洋观测装备主要包括科考船、浮标以及卫星等，但由于在海况适应性、可控性、电磁信号传输等方面存在局限性，难以针对具体台风核心区上层水体精准观测^[10]。随着自动化和信息技术的发展，近年来使用海洋机器人进行海洋动态精准观测的频次显著提高^[11]。但由于台风区域内天气恶劣、海况复杂，传统海洋机器人均存在高海况适应性差问题，无法获得有效观测数据。因此，亟需新的观测技术手段，以获取台风海气界面和混合层水体现场观测资料。

基于放风筝运动启发，结合海洋机器人技术和船舶拖曳技术提出一种由水面智能拖体、水下高机动AUV和高强度拖缆组成的“水下风筝”式三体AUV平台观测方案，如图1所示。拖体类似“水下风筝”式AUV中的“风筝”，主体外形为鱼雷形，两侧对称布置由电机驱动的翼板，可根据需求主动控制调节翼面角度以改变拖体水动力特性，拖体上通过搭载海表飞沫测量仪，可对海表飞沫的直径、数量等关键参数进行测量；拖缆是AUV和拖体之间的“牵引线”，通过分布式集成温度链实现对上层水体温度的测量；AUV是整个系统的动力源，通过拖缆牵引拖体运动，同时可通过内置自动绞车调节缆长，结合拖体的姿态调节功能协同控制拖体在水面跟随波浪起伏前向运动而不完全浸没水中（保持半潜状态），如图2所示。通过以上创新性设计，“水下风筝”式AUV有望实现台风核心区海表飞沫和上层水体温度的同步测量。

但台风区域海况恶劣，通常伴随着大风大浪等强天气因素干扰，“水下风筝”式AUV在极端动态海洋环境下的安全航行和观测功能可行性问题需要探究，因此将针对当前AUV拖曳系统的相关研究工作进行调研综述，主要包括AUV拖曳磁力计、AUV拖曳水听器阵列和AUV拖曳通信浮标3种。

AUV拖曳磁力计：由于船舶拖曳磁力计成本高、受天气条件限制且设备安全性差，使用AUV拖曳磁力计进行地磁测量成为AUV拖曳系统的应用领域之一。Dhanak等^[12]通过AUV拖曳10 m长拖缆和磁力计获得了较好的地磁测量数据，但航行时由于水动力作用导致磁力计、拖缆与AUV几乎处于同一水平面，无法应用于剖面水体观测，其高海况下运动性能也未得到验证。在此基础上Miller等^[13]利用Matlab-Orcaflex接口对AUV拖曳磁力计进行模拟，确定了最佳缆长和拖曳速度，但仍未达到剖面水体观测要求。Kostenko等^[14]分析拖曳式磁力计对AUV的扰动，建立了磁力计干扰数学模型，提出了控制补偿算法并通过初步实验验证了有效性，但其仍然存在缆长过短以及未考虑波浪影响问题。

AUV拖曳水听器阵列：欧盟“H2020计划”提出可扩展的广泛移动水下声呐技术（WiMUST），旨在通过协作海洋机器人系统实现用于地震勘测的分布式拖曳声学阵列技术^{[15 - 16]}。WiMUST其中很重要的一项工作即是使用AUV拖曳水听器阵列进行物探活动 ^[17]。WiMUST中使用2台AUV，每台拖曳8 m长水听器拖缆，分别在3 m和5 m深度航行，并成功获得大面积地震采集数据，但并未对动态海洋环境下AUV拖曳水听器阵列的运动进行分析^[18]。

AUV拖曳通信浮标：拖曳式浮标与AUV组成的AUV拖曳浮标系统在海底环境观测方面应用广泛，浮标用于通信和定位，通过对柔性拖缆方程求解，可得AUV相对浮标位置，以此进行AUV的定位与控制^[19]。Bykanova等^[20]设计了一种受到较小水动力阻力作用的AUV拖曳式浮标，并通过实验验证了其可行性。但实验中AUV下潜深度不超过10 m，拖曳缆长不超过20 m，且并未考虑波浪对浮标、拖缆及AUV的影响，不适用于动态海洋剖面观测。Pan等^[21]基于集中质量法、力矩定理和角动量定理，建立了浮标和拖缆的动力学模型，模拟了AUV发射拖曳浮标过程，得到能够较好预测浮标轨迹和拖缆形状的数学模型，该模型同样未考虑波浪对浮标及拖缆影响。

因此国内外针对AUV拖曳系统的研究多集中于静水情况下的力学分析，而高海况条件下，相对于传统AUV拖曳系统，“水下风筝”式AUV存在AUV力学响应不清。因此本文以“水下风筝”式AUV在波浪下的力学响应特性为研究目标，首先建立基于微振幅波理论和经典多体刚柔结构动力学建模方法的“水下风筝”式AUV波浪影响下的动力学模型，仿真模拟了规则波不同位置、不同拖曳速度、不同海况等级对AUV端拉力、深度和所需缆长的影响，为从当前海洋机器人技术现状评估台风条件下“水下风筝”式AUV海气界面与上层水体同步观测方案可行性提供重要判定依据。

1 动力学建模 1.1 坐标系建立与受力分析

为便于研究，本文仅对“水下风筝”式AUV在垂直面内受波浪影响的力学特性进行分析。定义惯性坐标系$ {O}-{X}{Y} $、拖缆微元局部坐标系$ {{O}}_{1}-{{X}}_{1}{{Y}}_{1} $如图3所示，图中$ {O}{X} $轴与$ {{O}}_{1}{{X}}_{1} $轴夹角为$ \theta $，U为AUV航行速度，U_t为拖体速度，H为有义波高，$ \lambda $为波长，h为AUV相对静水深度，所用单位均为国际标准单位。

“水下风筝”式AUV受力情况如图4所示。图中$ \boldsymbol{i} $、$ \boldsymbol{k} $为惯性坐标系单位矢量，$ \boldsymbol{t} $、$ \boldsymbol{n} $为拖缆坐标系单位矢量。拖缆与AUV连接处为首端，与拖体连接处为末端，$ \boldsymbol{T} $为拖缆张力，$ {\boldsymbol{T}}_{0} $、$ {\boldsymbol{T}}_{\boldsymbol{n}} $分别为AUV和拖体受到来自拖缆的拉力。惯性坐标系单位向量和局部坐标系单位向量的变换为：

$ \left\{\begin{array}{c}\boldsymbol{t}=\boldsymbol{i}\mathrm{c}\mathrm{o}\mathrm{s}\theta +\boldsymbol{k}\mathrm{s}\mathrm{i}\mathrm{n}\theta，\\ \boldsymbol{n}=\boldsymbol{k}\mathrm{c}\mathrm{o}\mathrm{s}\theta -\boldsymbol{i}\mathrm{s}\mathrm{i}\mathrm{n}\theta。\end{array}\right. $

(1)

1.2 微振幅波理论

微振幅波理论是研究各种海洋问题时必不可少且最为基本的波浪理论。根据微振幅波理论，波面以简谐形式进行起伏运动，波浪内的水质点以封闭椭圆为运动轨迹作圆周运动，波浪中线与静水面重合。水面波动可以描述为：

$ Y=\frac{H}{2}\mathrm{s}\mathrm{i}\mathrm{n}\left(kX-\omega t\right)。$

(2)

式中：$ H $为有义波高；$ \omega $为圆频率；$ k $为波数；$ k=2{\text{π}} /\lambda $，$ \lambda $为波长，$ {h}_{0} $为海底与静水面之间的垂直距离。

当$ {h}_{0}\geqslant \lambda /2 $时，产生的波浪称为深水波，当$ {h}_{0}\to \infty $时，可得深水波速度场表达式^[22]:

$ \left\{\begin{aligned} & {J}_{X}=\frac{kHg}{2\omega }{e}^{kY}\mathrm{c}\mathrm{o}\mathrm{s}(kX-\omega t)，\\ & {J}_{Y}=\frac{kHg}{2\omega }{e}^{kY}\mathrm{s}\mathrm{i}\mathrm{n}(kX-\omega t)。\end{aligned}\right. $

(3)

式中：$ {J}_{X} $、$ {J}_{Y} $分别为波浪内不同位置水质点运动速度沿惯性坐标系$ {O}{X} $、$ {O}{Y} $轴的分量，$ X $、$ Y $为拖缆微元在惯性坐标系中的位置，$ X'=\left(1+eT\right)\cos\theta $，$ Y'=(1+ eT)\sin\theta $。

由$ {\boldsymbol{J}=J}_{X}\boldsymbol{i}+{J}_{Y}\boldsymbol{k}={J}_{t}\boldsymbol{t}+{J}_{n}\boldsymbol{n} $可得：

$ \left\{\begin{array}{c}{J}_{t}={J}_{X}\mathrm{cos}\left(\theta \right)+{J}_{Y}\mathrm{sin}\left(\theta \right)，\\ {J}_{Y}={J}_{Y}\mathrm{cos}\left(\theta \right)-{J}_{X}\mathrm{sin}\left(\theta \right)。\end{array}\right. $

(4)

式中，$ {J}_{t} $、$ {J}_{n} $为波浪内不同位置水质点运动速度沿局部坐标系$ {{O}}_{1}{{X}}_{1} $、$ {{O}}_{1}{{Y}}_{1} $轴的分量。

1.3 动力学方程

拖体与拖缆末端相连，可通过改变翼面角度改变自身升阻特性。由图4分析可知，拖体受到的力分别有拖缆拉力$ {\boldsymbol{T}}_{\boldsymbol{n}} $，流体升力$ {\boldsymbol{F}}_{\boldsymbol{L}} $、流体阻力$ {\boldsymbol{F}}_{\boldsymbol{D}} $、重力$ {\boldsymbol{F}}_{\boldsymbol{g}} $、浮力$ {\boldsymbol{F}}_{\boldsymbol{v}} $、$ \phi $为波面角，$ \theta $为拖缆末端与水平方向的夹角。由于拖体随着波浪起伏，且尺寸远小于波长，为便于研究，仅计算一阶波浪力，并将拖体视作紧贴波面运动的质点。

波浪中拖体受力情况可表示为：

$ \boldsymbol{T}_{\boldsymbol{n}}+\boldsymbol{F}_{\boldsymbol{D}}+\boldsymbol{F}_{\boldsymbol{L}}+\boldsymbol{F}_{\boldsymbol{g}}+\boldsymbol{F}_{\boldsymbol{v}}=m_t\boldsymbol{a}_{\boldsymbol{t}}。$

(5)

式中：$ {m}_{t} $为拖体质量；$ {\boldsymbol{a}}_{\boldsymbol{t}} $为拖体加速度。拖体在波浪中航行时受到的流体力主要为升力$ {\boldsymbol{F}}_{\boldsymbol{L}} $和阻力$ {\boldsymbol{F}}_{\boldsymbol{D}} $，其表达式为：

$ \left\{\begin{aligned} & {F}_{L}={\frac{1}{2}C}_{L}{\rho {U}_{t}}^{2}{S}_{t}，\\ & {F}_{D}={\frac{1}{2}C}_{D}{\rho {U}_{t}}^{2}{S}_{t}。\end{aligned}\right. $

(6)

式中：$ \rho $为海水密度；$ {S}_{t} $为拖体迎流面积；$ {C}_{L} $、$ {C}_{D} $分别为拖体的升力系数和阻力系数。

将拖体受力在$ \boldsymbol{i} $、$ \boldsymbol{k} $方向上分解可得：

$ \boldsymbol{T}\mathrm{c}\mathrm{o}\mathrm{s}\left(-\theta\right)-\boldsymbol{F_L}\mathrm{s}\mathrm{i}\mathrm{n}\phi-\boldsymbol{F_D}\mathrm{c}\mathrm{o}\mathrm{s}\phi=0，$

(7)

$ \boldsymbol{F_v}+\boldsymbol{F_L}\mathrm{c}\mathrm{o}\mathrm{s}\phi-\boldsymbol{F_g}-\boldsymbol{T}\mathrm{s}\mathrm{i}\mathrm{n}\left(-\theta\right)-\boldsymbol{F_D}\mathrm{s}\mathrm{i}\mathrm{n}\phi=m_ta_t。$

(8)

当拖体处于波浪中不同位置时，由以上方程可得拖缆拉力大小$ \boldsymbol{T_n} $和拉力角度$ \theta $。

拖缆首端连接至AUV，末端连接至拖体，由AUV拖曳前进，前进过程中拖缆受到重力、浮力、水动力和惯性力共同作用。拖缆上任意微元的受力情况可表示为：

$ \frac{\partial }{\partial S}\boldsymbol{T}+\boldsymbol{W}+\boldsymbol{F}+\boldsymbol{B}=0。$

(9)

$ \boldsymbol{T} $为拖缆微元上的张力，方向沿拖缆切向：

$ \boldsymbol{T}=T\boldsymbol{t} 。$

(10)

$ \boldsymbol{W} $为拖缆微元受力拉伸后的负浮力，方向垂直向下：

$ \boldsymbol{W}=-\left(m-\rho A\right)g\boldsymbol{k}/{S}'=-{w}_{1}\boldsymbol{k}/S'。$

(11)

式中：$ m $为未拉伸时拖缆单位长度的质量；$ \rho $为海水密度；$ A $为未拉伸时拖缆的横截面积；$ {w}_{1} $为未拉伸时单位长度拖缆在海水中的重量，$ {w}_{1}=\left(m-\rho A\right)g $；$ g $为重力加速度；$ S $为拖缆微元受力拉伸后的长度，$ S=(1+ \varepsilon)s $，$ s $为拖缆微元未拉伸时的长度，$ \varepsilon $为应变，$ \varepsilon=eT $、$ e=1/EA $，$ E $为拖缆的杨氏模量。

将$ \boldsymbol{W} $在$ \boldsymbol{t} $、$ \boldsymbol{n} $两方向分解得到：

$ \boldsymbol{W}=-{w}_{1}(\boldsymbol{t}\mathrm{s}\mathrm{i}\mathrm{n}\theta +\boldsymbol{n}\mathrm{c}\mathrm{o}\mathrm{s}\theta )/S'。$

(12)

$ \boldsymbol{F} $为单位长度拖缆受力拉伸后在水中受到的阻尼力，对截面为圆形的拖缆，阻尼力表示为^[23]：

$ \begin{split} \boldsymbol{F}=\,& -\frac{1}{2}\rho d\left\{{\text{π}} {C}_{t}\left(\boldsymbol{U}\cdot \boldsymbol{t}\right)\left|\boldsymbol{U}\cdot \boldsymbol{t}\right|\boldsymbol{t}+ \right.\\ & \left.{C}_{n}\left|\boldsymbol{U}-\left(\boldsymbol{U}\cdot \boldsymbol{t}\right)\boldsymbol{t}\right|\left[\boldsymbol{U}-\left(\boldsymbol{U}\cdot \boldsymbol{t}\right)\boldsymbol{t}\right]\right\}。\end{split} $

(13)

式中：$ d $为受力拉伸后的拖缆直径；$ {C}_{t} $为拖缆切向阻力系数；$ {C}_{n} $为拖缆法向阻力系数。

将$ \boldsymbol{F} $在$ \boldsymbol{t} $、$ \boldsymbol{n} $两方向分解得到：

$ \boldsymbol{F}=-\frac{1}{2}\rho d\left({\text{π}} {C}_{t}{V}_{t}\left|{V}_{t}\right|\boldsymbol{t}+{C}_{n}{V}_{n}\left|{V}_{n}\right|\boldsymbol{n}\right) 。$

(14)

$ \boldsymbol{B} $为受力拉伸后单位长度拖缆在水中受到的惯性力，包括达朗贝尔惯性力和理想流体惯性力。

$ \boldsymbol{B}=-\frac{\partial }{\partial t}\left\{\left[{m}_{1}\dot{\boldsymbol{r}}-\rho A\boldsymbol{J}-\rho A\left(\boldsymbol{U}\cdot \boldsymbol{t}\right)\boldsymbol{t}\right]/S'\right\}。$

(15)

式中，$ \boldsymbol{J} $为波浪在拖缆微元处的水质点速度，由于$ \boldsymbol{J}={\boldsymbol{J}}_{\boldsymbol{i}}+{\boldsymbol{J}}_{\boldsymbol{k}}={\boldsymbol{J}}_{\boldsymbol{t}}+{\boldsymbol{J}}_{\boldsymbol{n}} $，将$ \boldsymbol{B} $在$ \boldsymbol{t} $、$ \boldsymbol{n} $两方向分解得到：

$ \boldsymbol{B}=-\frac{\partial }{\partial t}\left\{\left[m{V}_{t}\boldsymbol{t}+{m}_{1}{V}_{n}\boldsymbol{n}-\rho A{J}_{t}\boldsymbol{t}-\rho A{J}_{n}\boldsymbol{n}\right]/S'\right\} $

(16)

进一步计算得到：

$ \begin{split} & {\boldsymbol{B}=}\\ & {-\left(\frac{m\dot{{V}_{t}}\boldsymbol{t}+m{V}_{t}\dot{\theta }\boldsymbol{n}+{m}_{1}\dot{{V}_{n}}\boldsymbol{n}-{m}_{1}{V}_{n}\dot{\theta }\boldsymbol{t}-\rho A\dot{{J}_{t}}\boldsymbol{t}-\rho A{J}_{t}\dot{\theta }\boldsymbol{n}-\rho A\dot{{J}_{n}}\boldsymbol{n}+\rho A{J}_{n}\dot{\theta }\boldsymbol{t}}{1+eT}\right)+}\\ & {\frac{\boldsymbol{e}\dot{\boldsymbol{T}}\left(\boldsymbol{m}{\boldsymbol{V}}_{\boldsymbol{t}}\boldsymbol{t}+{\boldsymbol{m}}_{1}{\boldsymbol{V}}_{\boldsymbol{n}}\boldsymbol{n}-\boldsymbol{\rho }\boldsymbol{A}{\boldsymbol{J}}_{\boldsymbol{t}}\boldsymbol{t}-\boldsymbol{\rho }\boldsymbol{A}{\boldsymbol{J}}_{\boldsymbol{n}}\boldsymbol{n}\right)}{{\left(1+\boldsymbol{e}\boldsymbol{T}\right)}^{2}}}。\\[-1pt] \end{split} $

(17)

式中：$ \boldsymbol{U} $为拖缆微元相对海水的速度，$ \boldsymbol{U}={V}_{t}\boldsymbol{t}+{V}_{n}\boldsymbol{n} $；$ \boldsymbol{r} $为拖缆微元的位置矢量。

在“水下风筝”式AUV中，局部坐标系$ {{O}}_{1}-{{X}}_{1}{{Y}}_{1} $的单位矢量$ \boldsymbol{t} $、$ \boldsymbol{n} $的偏微分可以表示为^[24]：$ \dot{\boldsymbol{t}}=\dot{\theta }\boldsymbol{n} $，$ \dot{\boldsymbol{n}}=-\dot{\theta }\boldsymbol{t} $，$ \boldsymbol{t}'=\theta '\boldsymbol{n} $，$ {\boldsymbol{n}}'=-\theta '\boldsymbol{t} $。

根据拖缆连续性方程：

$ \left(\boldsymbol{r}'\right)\dot{}=\left(\dot{\boldsymbol{r}}\right)' $

(18)

$ \frac{\partial }{\partial \boldsymbol{t}}\left[\left(1+\boldsymbol{e}\boldsymbol{T}\right)\boldsymbol{t}\right]=\frac{\partial }{\partial \boldsymbol{s}}\left({\boldsymbol{V}}_{\boldsymbol{t}}\boldsymbol{t}+{\boldsymbol{V}}_{\boldsymbol{n}}\boldsymbol{n}\right) $

(19)

$ \boldsymbol{e}\dot{\boldsymbol{T}}\boldsymbol{t}+\left(1+\boldsymbol{e}\boldsymbol{T}\right)\dot{\boldsymbol{\theta }}\boldsymbol{n}={\boldsymbol{V}}_{\boldsymbol{t}}'\boldsymbol{t}+{\boldsymbol{V}}_{\boldsymbol{n}}'\boldsymbol{n}+{\boldsymbol{V}}_{\boldsymbol{t}}\boldsymbol{\theta }'\boldsymbol{n}-{\boldsymbol{V}}_{\boldsymbol{n}}\boldsymbol{\theta }'\boldsymbol{t}。$

(20)

式（17）、式（20）在$ \boldsymbol{t} $、$ \boldsymbol{n} $方向进行分解，联立式（7）～式（8）以及拖缆微元在惯性坐标系中的位置，可得到受波浪影响的拖缆动力学偏微分方程组，写成以下矩阵形式：

$ \boldsymbol{M}{\boldsymbol{Y}}'=\boldsymbol{N}\dot{\boldsymbol{Y}}+\boldsymbol{Q} 。$

(21)

式中：$ \boldsymbol{Y}={\left[\begin{array}{cccccccc}T\; \; {V}_{t}\; \; {V}_{n}\; \; \theta \; \; {J}_{t}\; \; {J}_{n}\; \; X\; \; Y\end{array}\right]}^{{\mathrm{T}}} $；$ {\boldsymbol{Y}}'=\dfrac{\partial Y}{\partial s} $；$ \dot{Y}=\dfrac{\partial Y}{\partial t} $。

$ \boldsymbol{M}=\left[\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& {-V}_{n}& 0& 0& 0& 0\\ 0& 0& 1& {V}_{t}& 0& 0& 0& 0\\ 0& 0& 0& T& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]，$

$ \scriptsize\boldsymbol{N}=\left[\begin{array}{cccccccc}\frac{-e\left(m{V}_{t}-\rho A{J}_{t}\right)}{1+eT}& m& 0& {m}_{1}{V}_{n}-\rho A{J}_{n}& -\rho A& 0& 0& 0\\ e& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1+eT& 0& 0& 0& 0\\ \frac{-e\left({m}_{1}{V}_{n}-\rho A{J}_{n}\right)}{1+eT}& 0& {m}_{1}& m{V}_{t}-\rho A{J}_{t}& 0& -\rho A& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right] ，$

$ \boldsymbol{Q} = \left[ \begin{array}{c}{w}_{1}\mathrm{s}\mathrm{i}\mathrm{n}\theta +\dfrac{1}{2}\rho D{\text{π}} \sqrt{1+eT}{C}_{t}{V}_{t}\left|{V}_{t}\right|\\ 0\\ 0\\ {w}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\theta +\dfrac{1}{2}\rho D\sqrt{1+eT}{C}_{n}{V}_{n}\left|{V}_{n}\right|\\ \dfrac{kHg}{2\omega }{e}^{kY}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left(kX-\omega t\right)\mathrm{c}\mathrm{o}\mathrm{s}\theta +\mathrm{s}\mathrm{i}\mathrm{n}\left(kX-\omega t\right)\mathrm{s}\mathrm{i}\mathrm{n}\theta \right]-{J}_{t}\\ \dfrac{kHg}{2\omega }{e}^{kY}\left[\mathrm{s}\mathrm{i}\mathrm{n}\left(kX-\omega t\right)\mathrm{c}\mathrm{o}\mathrm{s}\theta -\mathrm{c}\mathrm{o}\mathrm{s}\left(kX-\omega t\right)\mathrm{s}\mathrm{i}\mathrm{n}\theta \right]-{J}_{n}\\ \left(1+eT\right)\mathrm{c}\mathrm{o}\mathrm{s}\theta \\ \left(1+eT\right)\mathrm{s}\mathrm{i}\mathrm{n}\theta \end{array} \right] 。$

拖缆首端与AUV相连，且在AUV上连接位置固定，因此首端速度与AUV拖曳速度相同，且$ \boldsymbol{U}={V}_{t}\boldsymbol{t}+{V}_{n}, \boldsymbol{n}={V}_{X}\boldsymbol{i}+{V}_{Y}\boldsymbol{k} $，故首端边界条件可表示为：

$ \left\{\begin{array}{c}{V}_{t}={V}_{X}\mathrm{c}\mathrm{o}\mathrm{s}\theta +{V}_{Y}\mathrm{s}\mathrm{i}\mathrm{n}\theta，\\ {V}_{n}={V}_{Y}\mathrm{c}\mathrm{o}\mathrm{s}\theta -{V}_{X}s\mathrm{i}\mathrm{n}\theta 。\end{array}\right. $

(22)

拖缆末端与拖体相连，由式（7）～式（8）可得拖缆拉力$ {T}_{n} $和拉力角度$ \theta $，以此作为拖缆动力学偏微分方程组末端边界条件。

2 仿真计算结果及分析

本文使用Matlab和Comsol Multiphysics软件系数型偏微分方程组模块，采用有限差分法对建立的方程组进行时间和空间离散，联立首端边界条件和末端边界条件，使非线性偏微分方程组封闭，采用牛顿迭代法进行数值仿真，分析了规则波中不同位置、不同海况等级、不同拖曳速度下AUV端拉力和深度的变化情况。

本文所研究的“水下风筝”式AUV设计方案中拖体物理模型和相关参数分别如图5和表1所示。

拖缆相关参数如表2所示，由于拖缆在垂直平面内计算，故不考虑拖缆扭转变形。

不同海况时波浪参数如表3所示。

2.1 拖缆姿态及张力沿拖缆变化

在“水下风筝”式AUV执行高海况下的剖面观测任务时，拖缆需张紧且穿过上层水体，且观测过程中不发生断裂。因此拖缆姿态及沿拖缆张力的计算十分重要。

仿真条件设定为8级海况，AUV相对静水面定深100 m且分别以2 kn、3 kn和4 kn的速度前进。由于拖体随波面起伏，选取一个周期内拖体处于波峰、波谷和波浪中线4个位置进行仿真。根据拖体所处位置不同，得到拖缆姿态如图6所示，沿拖缆的张力变化情况如图7所示。

由图6可知，当AUV以不同速度前进时，拖缆在流体力和拖体产生的升力作用下，在水中均呈抛物线姿态。随着拖曳速度的增加，拖体与AUV的水平距离逐渐增大，当拖体位于波峰位置时水平距离最大，如图6（a）所示，2 kn、3 kn和4 kn时水平距离分别为743 m、776 m和787 m；拖体处于波浪中线下降位置时水平距离最小，2 kn、3 kn和4 kn时分别为315 m、362 m和393 m。当速度不变，由于拖体在波浪中位置不同，拖缆末端拉力和缆长均发生变化，引起拖体与AUV水平距离改变，当拖体从波浪中线到达波峰或波谷时，水平距离增加，当拖体从波峰或波谷到达波浪中线时，水平距离减少。

由图7可知，拖缆内部张力沿拖缆增加，在拖缆首端达到最大，其增大速率随拖曳速度的增加而增大。当AUV以不同速度航行，拖缆内各点张力随拖曳速度的增加而显著增大。当拖曳速度为2 kn时，拖缆内部张力变化范围为100 ~ 200 N，3 kn时内部张力变化范围为300 ~ 500 N，4 kn时拖缆内部张力变化范围为600 ~ 800 N，均未超过所选拖缆额定张力。

拖体位于波浪中不同位置时，拖缆内部张力不同，其中当拖体位于波浪中线时，拖缆内部张力最大，以图7（d）为例，拖曳速度分别为2 kn、3 kn和4 kn时，内部张力分别为195~218N、435~486N和757~851N。当拖体从波浪中线到达波峰或波谷时，沿拖缆张力整体减少，当拖体从波峰或波谷到达波浪中线时，沿拖缆张力整体增大。

结合图6和图7可知，8级海况时，在本文给定的任一拖曳速度下，“水下风筝”式AUV中的拖缆均能张紧且不发生断裂。随着拖体在波浪中位置改变，拖体受到的流体升力与阻力发生变化，引起拖缆作用在拖体上的拉力和拉力角度变化，导致拖缆姿态及沿拖缆张力发生改变。

2.2 拖曳速度影响分析

为进一步研究不同拖曳速度对AUV端受到拉力的影响，在8级海况下，分别使AUV深度和缆长固定，研究分析AUV端拉力变化情况。为使结果更具可信度，在一个波长内等间距取8个位置点，对拖体位于波浪中各位置时缆长、AUV深度及所受拉力进行仿真计算

当AUV相对静水面深度为100 m时。仿真结果如图8所示，其中${X} $轴为拖体位于一个波长内的位置，由于不同海况波长不同，${X} $轴以拖体处于波浪中线上升位置为零点，${X}=1\mathrm{\lambda } $时即为一个波长后同一位置。

由图8（a）可知，当海况等级不变，拖曳速度增大时，维持AUV航行深度所需的拖缆长度增加。当拖体处于波谷位置，不同拖曳速度所对应的缆长近乎相等。当拖体从波浪中线到达波峰或波谷时，缆长增加，当拖体从波峰或波谷到达波浪中线时，缆长减小。

由图8（b）可知，当拖体处于波浪中任一位置时，AUV端拉力受到的拉力均随拖曳速度的增加而增大，2 kn时拉力范围为162.45 ~ 223.28 N，3 kn时为408.69 ~ 497.65 N，4 kn时为752.2 ~ 872.67 N。当拖曳速度为2 kn、3 kn和4 kn时，AUV端拉力差分别为60.83 N、88.96 N和120.47 N，占最小拉力的37%、21.7%和16%。

由于AUV端受到拉力大小取决于拖缆末端张力，且拖体受到的流体升力和阻力与拖体速度的平方正相关，根据式（5），可知随拖曳速度增加，拖缆末端张力增大，故AUV端拉力显著增大。同理，当拖曳速度固定，拖体从波浪中线到达波峰或波谷时，其相对海水的速度减小，受到的流体升力和阻力降低，拖缆末端张力降低，故AUV端受到拉力变小。

当缆长为800 m时，AUV深度和受到拉力仿真结果如图9所示。

可知，拖缆长度为800 m时，AUV深度能够满足观测要求，拖曳速度增大时，AUV受到的拉力增大，且其相对静水面深度减少。拉力变化趋势与AUV深度100 m时的仿真结果相同。当拖体从波浪中线到达波峰或波谷时，其受到的流体力减小，拖缆末端张力降低，AUV端受到拉力降低且AUV深度减小。在此基础上，当拖体处于波峰位置时，拖缆末端高出静水面，此时AUV深度最小。由此可知，当缆长不变时，为了满足上层水体观测需求，AUV需要根据拖体处于波浪中的位置来调整深度，以确保拖体能够维持海面的运动。

2.3 AUV深度影响分析

考虑当前海洋机器人技术条件，拖曳速度定为2 kn，在8级海况下，仿真计算AUV相对静水面深度为80~100 m时拖缆和AUV端受到拉力，结果如图10所示。

可知，当AUV相对静水面深度增加，缆长和AUV端拉力均增加。当AUV深度分别为80 m、90 m和100 m，拖体处于波浪中线上升位置时缆长最小，分别为203 m、269 m和350 m，拖体处于波峰位置时缆长最大，分别为532 m、648 m和763 m。由前文仿真极端结果可知，拖缆呈抛物线姿态，故AUV深度增加时，缆长增加较为明显。

由于拖曳速度固定，拖缆末端张力和角度不变，AUV端受到拉力变化仅和缆长有关且变化程度较小，3种AUV深度下，根据拖体所处位置不同，AUV受到拉力均在140 ~ 230 N之间变化，符合现有海洋机器人承载能力范围。

2.4 海况等级影响分析

AUV相对静水面深度为100 m，拖曳速度为2 kn，仿真计算不同海况对缆长及AUV端拉力的影响，结果如图11所示。

可知，当海况等级分别为4级、6级和8级时，拖缆长度变化范围为490 ~ 719 m、412 ~ 737 m和350 ~ 763 m。当拖体处于波峰附近位置时，缆长达到最大，且缆长随着海况等级的增加而增大。除波峰附近位置外，拖缆长度随着海况等级增加而减少，当拖体处于波浪中线上升位置时，缆长最小，4级、6级和8级海况时分别为490 m、412 m和350 m。海况等级提高时，有义波高增加，当拖体位于波峰附近位置时，AUV与拖体垂直距离大于AUV相对静水面深度，因此所需缆长增大。

如图11（b）所示，当拖体处于波谷附近位置，AUV端拉力随着海况等级的增加而降低。除波谷附近位置外，拉力随着海况等级的增加而增加。当海况等级分别为4级、6级和8级时，AUV端拉力变化范围分别为153.14 ~ 189.59 N、155.37 ~ 205.66 N和158.63 ~223.28 N。

当拖体处于波浪中任一位置时，AUV端拉力差随海况等级增加而增大，8级海况时拉力差为64.65 N，占最小拉力的40.7%。

3 结　语

从本文研究结果来看，AUV拖曳系统涉及的静水或者单域（纯水域，不涉及跨海气界面）环境下的动力学问题国际上已经开展了一些研究，但针对高海况下的跨域运动响应问题还未有效解决。本文以“水下风筝”式AUV为研究对象结合微振幅波理论和经典多体刚柔结构动力学建模方法，建立了波浪影响下的“水下风筝”式AUV动力学模型，并以规则波中不同位置、不同海况等级、不同拖曳速度为自变量，仿真分析了AUV端所受拉力大小和所处深度值。从分析结果来看，“水下风筝”式AUV在高海况下航行具有可行性，将为台风核心区海气界面和上层水体同步动态观测提供可能，具有较大潜在应用价值。