下放式海洋微结构湍流剖面仪运动学与动力学分析*

2019-10-12宋大雷王浩杰周丽芹臧胜波

中国海洋大学学报(自然科学版) 2019年11期

宋大雷,王浩杰,周丽芹,臧胜波

(中国海洋大学工程学院，山东青岛 266100)

海洋湍流剪切数据是进行海洋微结构湍流研究的重要原始资料，对建立海洋物理模型进而研究海洋宏观运动有重要意义[1]。通过剪切数据计算得到的湍流动能耗散率，是湍流的一个重要特征参数。

Osborn团队于1970年代首先研制出搭载翼型剪切流传感器的垂直微结构湍流剖面仪，并成功进行测量海洋微结构湍流[2]。垂直微结构湍流剖面仪(以下简称湍流仪)是目前获取海水剪切流数据最普遍和最有效的平台之一，已受到海洋领域科研人员的普遍认可。

加拿大Rockland公司研制的VMP系列和德国ISW Wasser公司与SST公司合作研发MSS系列湍流仪是当前最具代表性的，得到国际市场的广泛认可。在国内，天津大学率先开展对湍流仪的研制，目前已研制成功第三代样机。本文以实验室自主研发的下放式湍流仪为研究对象，对其相关性能进行分析。

湍流仪下潜过程中，剪切传感器通过探针以一定的频率测量微结构湍流水平方向的脉动信号。根据下潜速度和攻角，计算微结构湍流脉动在时域、频域和空间波数域的强度，进而计算湍流动能耗散率[3]。湍流脉动频率和下潜速度成正比，而剪切传感器的固有频率决定脉动频率的测量上限[4]。剪切信号在频域和波数域上的分析皆基于湍流的脉动频率。

姿态会影响水流攻角，在小范围内，剪切传感器输出信号幅值与水流攻角近似成正比[5]。过大的倾角会导致测量误差加剧，给数据分析带来困难。

仪器的固有频率决定流致振动的频率[6]，而流致振动频率决定湍流耗散率的测量下限[7]。同时，流致振动会带来严重的噪音信号，影响剪切传感器的工作性能。

因此，合理而稳定的下潜速度、良好的下潜姿态以及小固有频率，对于准确可靠的测量湍流动能耗散率具有十分关键的作用。通过动力学分析，研究结构对这些因素间影响，能够为后期数据处理及结构优化提供帮助。

1 运动模型建立

1.1 运动坐标系

湍流仪在下潜的时候，由于受到流场的的影响，受力不平衡，会出现绕轴摆动。此时，其速度可以分解为两部分：旋转中心的直线速度，以及绕轴转动。

为便于分析，参考水下滑翔机的运动模型建立方法[8]，采用双坐标系对湍流仪的运动进行分析，坐标系如图1所示。

图1 运动坐标系Fig.1 Motion coordinate systems

图1中，OXYZ是绝对坐标系，O为仪器入水点，XOY为海平面，Z为竖直向下方向。oxyz为相对坐标，以仪器自身为参考，o为仪器自由转动的中心点，x为轴线方向，xoy为XOZ的平行面。

湍流仪内置姿态传感器，用偏航角、俯仰角和翻滚角(Heading-Pitch-Roll)三个角参数来描述仪器姿态。在OXYZ坐标系下，描述仪器的三个姿态角，坐标系如图2所示。

图2 姿态角参考坐标系Fig.2 Coordinate system for attitude angles

在图2中，OP为湍流仪轴线的平行线，OP1和OP2分别为OP在YOZ′和XOZ平面上的投影。β、φ、γ分别对应偏航角、俯仰角和翻滚角[9]。由于MTI姿态传感器安装时固定面的不同，三个姿态角数据可能存在互换。

角θ为仪器的倾角，即与竖直方向的夹角。仪器在平面内摆动时，β、φ符号要么保持相同，要么不同，变化率也成线性关系。因此，定义θ与β同正负，其值可以通过角β、φ计算得到：

(1)

在倾斜较小的情况下，上式通过近似简化为：

(1*)

1.2 动力学模型

在局部海域，海流方向具有一致性，即V在一个确定的平面内，假设其在Y轴(z轴)上的速度分量为0。为简化分析，将湍流仪在空间的运动降到二维坐标进行分析，其运动模型如图3。

图3 动力学模型图示Fig.3 Graphical kinetic model

在图3中，V为湍流仪o点相对周围流体平均流动的相对运动速度，m/s；χ为V与Z轴夹角，rad。G和B分别为重力和浮力，N；L1和L2分别为重心和浮心到旋转中心的距离，m。Fv为流体对仪器y轴方向的冲击力合力，N；Fd为流体对仪器x轴方向的阻力合力，N；Mw为流体作用力对o点的合力矩，N·m。不考虑湍流脉动，水流平均速度对剪切传感器的攻角α=θ-χ。

在二维坐标系下，湍流仪模型有三个自由度，即三个变量：X、Z和θ。此模型下，湍流仪的速度可分解为：

(2)

(3)

(4)

基于上述模型，在两个坐标系下分别建立动力学微分方程组：

(5)

(5*)

式中：m为总质量，kg；J为等效转动惯量，kg·m2。

2 参数计算

湍流仪总长为1 800 mm，耐压舱长度810 mm，外径125 mm，浮体长780 mm，外径D=220 mm。以仪器最尾端为参考平面，其与仪器轴线的交点为坐标原点，建立坐标系，用于其他重要结构参数分析，其坐标系见图4。

图4 重浮力配置坐标系Fig.4 Coordinate system for configuration of gravity and buoyancy

在此参考系下，相关结构参数汇总见表1。

表1 结构动力学参数Table 1 Kinetic parameters of the structure

2.1 旋转中心与转动惯量

根据能量最低原理，湍流仪受到外力作用时，会绕转动惯量最小的轴旋转[10]。由于转动时有附加质量(排开水的质量)存在，转动中心o不在质心处，而是位于重浮心之间，到重心的距离为：

(6)

式中：m为湍流仪质量，mw为排水质量，L为重浮心距离。通过计算，得到L1=0.092 m，转动中心o位于距尾端0.722 m处。

仪器对旋转中心o的转动惯量为：

(7)

式中：J1为仪器转动惯量，J2为附加质量(即排开的水)转动惯量。将相关参数代入，可得J=16.57 kg·m2。

2.2 流体作用力

在动力学模型中，湍流仪上有三个与流场相关的作用力：阻力Fd、冲击力Fv以及冲击力矩Mw。

Fd力方向与Vx相反，其中阻力系数CD已经通过计算得到。阻尼刷丝长为200 mm，CD=2.945，计算公式如下：

(8)

(9)

式中：L′为最尾端到o点的距离，L为总长，D(l)为该处直径，m。若仪器不存在转动，仪器受到的冲击力合力为：

(10)

Fd作用于仪器一周，对o取矩为0，因此合力矩Mw只与Fv有关。在分析转动惯量时，已将附加质量考虑在内，此处只计算由Fv引起的力矩，其公式为：

(11)

式中，L″为最尾端到Vy=0点的距离，m。若仪器不存在转动，那么：

(12)

3 下潜速度

分析下潜速度不考虑仪器的倾斜，即仪器只有Z轴方向的速度。根据公式(5*(b))，结合坐标转换公式，可以得出关于下潜速度VZ的关系式：

(13)

上式是带二次项的微分方程，难以析解。但是，在足够长时间后，重力、浮力和水阻力达到平衡，加速度为0，速度趋于稳定，那么公式13可化为：

(14)

加入配重结构，可以改变仪器下潜速度，配重使用铅块密度为11.3 g/cm3，每100 g铅块在海水中净重为91 g。加入配重质量Nkg，净重(0.06+0.91N)kg。结合阻力Fd，此时仪器的稳定速度：

(14*)

式中，g为重力加速度，取9.81 m/s2。

4 稳定姿态

当仪器处于动态稳定，线加速度和角加速度均为0，即要满足∂Vx/∂t=0，∂Vy/∂t=0，∂2θ/∂t2=0，公式(5)代入参数后得：

(15)

此方程组的唯一解是：θ=0，Vy=0，Vx=0.105。转换为XOZ坐标系下的速度：

此即公式(14*)中，M=0时的下潜速度配置。此时，Z轴方向的速度，χ=0，α=0。换而言之，只要倾角不为0，仪器就不可能处于动态平衡，必然存在某一参数的加速度。

仅考虑姿态的平衡，∂2θ/∂t2=0，分析VX和VZ对θ的影响。将浮力用重力近似替代，根据公式(5(c))，代入力矩参数后：

(16)

(16*)

从上式可以得出，在VX和VZ确定的情况下，重浮心距离L越大，角θ越小。假定VZ=0，对θ和VX关系进行定性分析得：

(17)

其函数曲线见图5。

随着VX变大，θ几乎成二次型增大，VX对仪器稳定性影响十分明显。分析不同VZ下θ和VX和关系，具体结果见图6。

图5 VX对倾角的影响Fig.5 Effect to tilt angle from VX

图6 不同VZ对倾角的影响Fig.6 Effect to tilt angle from variational VZ

从结果可以看出，相同VX下，VZ越大，出现的倾角越小。因此，较大的下潜速度有助于仪器姿态的稳定。

倾角与攻角的差值χ只与VX和VZ有关：

(18)

相同VX下，VZ越大，χ越小，即姿态数据越接近于攻角。

综上可以得出，L和VZ越大，仪器稳定性越好，攻角也越小。

5 自由摆动

分析湍流仪下潜过程中的自由摆动，即θ的变化规律，假定V=0。从Mw可以看出，公式(5(c))是一个非线性振动模型的微分方程，Mw为与多个参数相关的非线性阻尼。若不考虑流场的作用力，只分析重力和浮力的作用。即不考虑阻尼的影响，取Mw=0。考虑在小范围摆动时(θ<5°)，sinθ近似于θ。那么可以将方程简化为线性无阻尼自由振动模型：

(19)

这是典型的简谐振动模型，结构固有频率为：

(20)

对于本结构，由于mw与m十分接近，用m代如计算。基于J1和J2不随L发生变化假定，上式可写为：

(20*)

其关系如7图所示。

图7 结构固有频率Fig.7 Natural frequency of structure in different configuration

从图7中可以看出，随着重浮心距离变大，固有频率也随之变大。将具体参数代入到公式中，可以得出：

(21)

在固有频率下，湍流仪倾角及角速度的时间函数为：

(22)

上式中，θm为摆动幅值，τ为初始相位。对应剪切传感器探针处因自由摆动带来y轴方向速度为：

(23)

在此摆动频率下，max(|at|) = 4.513θm，单位为m/s2。

无阻尼振动振幅不会衰减，实际中不存在。考虑阻尼的存在，根据公式11得：

(24)

因此完整的摆动模型为：

(25)

阻尼力大小与角速度平方成正比，是一个二次非线性阻尼。二次阻尼不会改变系统的固有频率，幅值衰减一段时间后系统开始做简谐振动[12]。通过数值算法来求近似解，定性分析振动衰减特性。定义t0时刻θ=0.087 rad(5°)，角速度为0。计算角度、角速度、角加速度随时间变化分别如图8所示：

图8 振动方程数值结算结果Fig.8 Numerical results of the vibration equation

图8给出了各个参数随时间的变化，可以看出，随着时间的推移，振动周期几乎不变，而振动幅值在阻尼力的作用下持续衰减。为了更明确的分析衰减特性，计算后统计前8个角度峰值出现的时间和幅值(见表2)。

表2 θ峰值数据Table 2 Data on peaks of tilt angle

衰减幅度η为与上一个幅值的比值。从表2中数据分析，系统的摆动周期为3.02 s，振动频率为0.331 Hz，即上面计算结构的固有频率。随着时间推移，幅值变动幅度越来越小，相应的阻尼系数也越来越小，直到振幅平稳。因此，在信号处理时，必须注意该频率振动的影响。

6 试验数据分析

6.1 配速试验

在湍流仪完成后，在青岛市南姜码头先后进行了多次配速试验。试验水深约20 m，相同配置进行多次试验，取平均值用于比较。摘取部分试验速度数据如表3所示。

表3 下潜速度试验数据摘录Table 3 Extract of test data on sinking velocity

将试验速度与理论速度进行比较，结果见图9。

图9 理论与试验结果比较Fig.9 Comparison of theoretical and experimental results

从图9中可以看出，在小于配置速度(<0.3 m/s)部分，试验速度明显低于理论速度，分析为仪器受缆绳拖曳影响较大所致。在配置速度范围(0.3～0.7 m/s)内，理论计算的速度误差较小，十分接近试验结果。之后，随着配置增大，试验速度变化率明显大于理论速度，可能是因为阻尼刷丝变形导致阻力系数变小。由于在码头试验受到风浪、堤坝和试验母船的影响，通过传感器数据计算以及试验时测得平均速度与实际速度间均可能存在较大误差。

6.2 姿态分析

2016年7月，湍流仪跟随东方红2号科考船在南海进行了第二次海试。海试中，将无包塑的细缆系于湍流仪尾部进行布放。将湍流仪下放至海面下后等待一段时间，待其稳定后松开缆绳同时让缆绳始终保持松弛状态。当仪器到达设定深度后，通过拖拽缆绳回收。此次海试共完成三个剖面的测量，成三个剖面的测量，相关信息记录如见表4。