深海潜标姿态的动力学仿真分析与验证及沉降原因探究

王 蓓 ,陈永华 ,于 非

(1.中国科学院海洋研究所,山东 青岛 266071; 2.中国科学院大学,北京 100049; 3.中国科学院海洋大科学研究中心,山东 青岛 266071)

现阶段,数字信息在水环境内很难有效传递[1],自容式海洋潜标布放到定点位置几个月至几年,待其回收之后才能获得布放时间段海洋资料数据[2-3]。20世纪 50年代,美国率先开展海洋潜标技术[4],带动了世界各个海洋强国潜标观测技术的发展,锚泊的设计和分析方法也随之进步[5]。

迄今为止,数学模型已经能够逐步模拟出极为近似锚泊线的柱状模型,不论是否采用弹性存在的方法,都能对其张力进行分析[6]。Smith等[7]采用拉格朗日迭代分析悬链线方程模拟深海锚系多种成分锚泊线,将问题转化为8自由度的单一多项式方程; Chai等[8]基于半解析准静态法模拟三维局部着陆和完全悬挂锚泊线的假想姿态; Xu等[9]应用ACER方法预测3小时极限系泊张紧力从而得出研究系泊缆可靠性的详细程序。计算方法的不断发展促使可视模型的广泛运用。Dewey[10]设计的Mooring design,可用于设计单点海洋系泊、评估风和流影响下的系泊张力和形状、模拟系泊受力时的元件位置; Shibata等[11]使用 OrcaF-lex计算电流发生器施加到动态电缆上的应力,通过将电缆固定在系泊绳的多个点位模拟出最安全样式。前者程序简单但预测精度低,模型单一; 后者多用于船舶工程中电路传送大系统的模拟。国内仿真方法各异,汪舟红等[12]为达成减少船舶运动、均衡缆绳张力的优化目标,运用 OPTIMOOR软件将不同的缆绳材质、长度以及预张力进行对比。李亚男等[13]运用水动力分析软件AQWA研究半潜式平台系统的定位方法,在风、浪、流联合作用下,通过改变缆长来模拟平台的定位。目前国内系泊仿真分析倾向于设备定位,缆系布放和其水下姿态对系泊影响虽多,研究较少。

深海潜标系统多为单点系泊[14],在无环境载荷时,其预张力使缆系呈铅直态,存在波和流的作用时,缆系发生不可避免的偏移沉降,作用力超过一定幅值会导致零张力和最大张力的出现[15]。本文所模拟的潜标布放于西太平洋黑潮延伸区[16],该位置具有独特的大、小双环流结构,区域内海流错综复杂,中尺度涡现象非常活跃[17],受其影响,海域流速比大洋的特征流速高数倍,且延伸至水下几百至上千米[18],大洋中的中尺度涡以与长Rossby波速相近的速度向西传播[19],其运输作用相应地改变了海水密度[20],促成小尺度能量的传递[21],是潜标环境载荷中不可忽略的一部分。故而为避免锚系受到异常海流的影响进而损坏,需对其设计提出相应的规避风险要求,以保证其在工作时得到的数据完整、真实、有效。为保证自容式潜标在水下工作状态的可预测性,同时能够最大程度的优化缆系的结构设计,以达到安全效应和经济效应相结合的目的,本文拟结合中国科学院海洋研究所在西太平洋(30°N,146°E)布放潜标的回收数据,通过深海三维环境的建立还原海流及中尺度涡作用的影响; 分析设备和缆绳的材料性质,整体建模还原深海锚泊潜标系统结构部件。参考辅助设计资料[22-24],应用ADAMS动力学仿真模型,将深海三维模型与锚系建模叠加,验证其分析系泊姿态的合理性以及在不同海流条件下的准确度,得出缆系张力结果,分析异常沉降形成原因,提高潜标布放安全性,达到观测所需要的稳定性。

1 深海锚泊系统动力学仿真原理

1.1 潜标系统

锚泊系统实际全长为5 870 m,布放水深6 070 m。该潜标系统主要包括1个主浮体、7组玻璃浮球组、2组温深盐链、深海海流观测单元、1组声学释放单元以及1组锚系单元,所有单元通过10 mm凯夫拉缆绳串联布放。主要搭载的设备如下: 主浮体搭载2台75 kHz声学多普勒剖面测量仪(ADCP),其流速观测范围为水面至水下1 000 m; 在温盐链的设计中,水深200 m至500 m均匀布放四台温盐深测量仪(SBE37CTD),水深250 m至1 000 m布放10台温仪(SBE56); 深海海流情况通过位于1 800 m、3 000 m和5 800 m的3个单点海流计观测深海的流速结构特征。潜标系统结构示意图如图1所示。

图1 缆系结构设置图Fig.1 Diagram of the cable system structure setting

1.2 潜标动力学分析方法

当整个锚系系统在有流海域布放时,水流冲击力、重力锚阻力、锚系重力和浮力会使整个锚系系统进行动态飘移,整个缆系必将在三维空间中沿两个水平的方向横摇 ΔX和 ΔY,在铅直方向上产生一定的沉降ΔZ[25],水流阻力方程如式(1)所示。

其中,ρ为流体密度,C1、C2分别为切向阻力系数和法向阻力系数,A1、A2为切向迎流面积和法向迎流面积,v为流动速度,θ为潜标锚泊设备与水平方向夹角。

缆体分析则是将柔性体绳分成有限个刚性有限元,使整个模拟系统形成以多体动力学理论为基础的集中质量-弹簧模型,其锚系缆绳上第i个质量的运动方程如式(2)所示[26]。

其中,mi为第i段缆绳的质量,ai为其加速度,ei+1/2和ei-1/2分别为结点i、i+1和i、i-1间被拖曳流体的虚质量,aiN|i+1/2和aiN|i-1/2是向量ai在两段上的法向分量,力向量Fi包括两段缆绳中的张力、拖曳力、重力和浮力,以及其他任何外力。

1.3 潜标锚系系统有限元分析

模拟仿真自容式潜标系统在水下工作时,将对系泊与环境进行如下假设,以便简化分析:

1) 缆绳模型由n个微分缆绳单元组成,缆绳单元为刚直杆体,两个微分单元之间进行铰接,传递受力。同时由于n足够大,铰接缆绳模型的净重以及张力弯曲性质等都无限接近实际缆绳模型。

2) 在全部深度上,海流方向都需分为u和v两个水平方向和一个垂直方向w,作用在任何一个元素上的三维模型方向一致,流体阻力作用在锚泊缆绳与水流速度向量的平面内。

3) 假定系泊连接缆绳为柔性结构,不能传递弯矩。

4) 永远假定锚系底部自由度为0。

假设中n可以无限大以达成精密度极高的模拟计算,但由于模型边界条件的理想化、回收数据量有限、忽略柔性体传递弯矩、水下震荡现象的作用等,可能导致模拟计算结果与实际情况存在一定差异。

根据假设绘制出第n个缆绳单元的受力分析图如图3所示,列出受力平衡姿态下计算公式。

图3 第n段锚泊缆绳受力图Fig.3 Force diagram of the nth section of the mooring rope

表1 锚系搭载设备水下净重力表Tab.1 Net weights of the underwater anchoring equipment

图2 集中质量-弹簧模型Fig.2 Lumped mass-spring model

式中,Tn+1和Tn-1分别为第n+1和第n-1个微元对于第n个微元产生的牵引力,这两个牵引力与铅直方向的夹角分别为θn+1和θn-1,将这两个牵引力在xz平面上的投影与x轴的夹角分别定为αn+1和αn-1,设B为微元n所受的浮力,G为微元所受的重力,而FU和FV分别为微元n在三维海流中受U向和V向海流影响所受到的黏性阻力与法向阻力之和,其中:

式中,Cn为缆绳法向/切向阻力系数,L为第n段微元长度,d为缆绳微元直径,vun、vvn分别为第n段微元所受U、V方向平均流速,θ为倾斜缆绳与平面xz之间夹角。

将FU、FV代入平衡姿态计算公式,由于第n个微元和第n+1个微元仍满足上述平衡方程,经过联立求解可知,当体系状态稳定,Tn-1、θn-1、αn-1以及海流情况vun、vvn确定时,第n段微元θ有且仅有唯一解,第n+1段微元Tn+1、θn+1、αn+1有且仅有唯一解。配合如下几何关系,通过分布外推法、打靶法等数值求解方法,即可得到锚泊线几何形状[27]。

1.4 潜标系统模型设计

据玻璃浮球和主浮体位置将锚系缆绳分成 8段进行等比例建模,由于不同深度存在不同流速,进而对于缆绳的精密度划分要求有明显区别,3 000 m以下海洋潜流影响极小的缆系受力单一且稳定,适当降低精密度可减少结构的复杂性。其中凯夫拉缆绳的参考水下密度为1.05 g/cm3。缆绳微元化分量n的取值如表2所示。

表2 缆绳微元划分表Tab.2 Cable micro-element divisions

缆绳建模完成后在结构相应位置固定等比例大小的CTD、单点海流计、ORE 等锚系悬挂设备,并施加同等重力。

2 潜标流阻环境载荷影响因素

2.1 潜标布放位置海流载荷影响

潜标锚系系统于2017年4月布放于日本东南部西太平洋(30°N,146°E),于 2018年 6月进行回收。潜标回收后获得 ADCP以及 SBE37CTD,提取设备升降情况和海洋要素,搭载双75 kHz ADCP的主浮体设备的深度变化情况,经过正交分解后的U、V和W三种方向时间-深度-海流关系图如图4、图5、图6。

图4 经向海流流速与主浮体深度图Fig.4 Diagram of the position relationship of the main floating body and zonal ocean-current velocity

图5 纬向海流流速主与浮体深度图Fig.5 Diagram of the position relationship of main floating body and meridional ocean current velocity

图6 垂向海流流速主与浮体深度图Fig.6 Diagram of the position relationship of the main floating body and vertical ocean current velocity

仅通过本次对垂向海流数据的处理可知布放期间海流方向铅直向下,依据公式(3)—公式(9)计算可知,在0～1 000 m内垂向海流对锚泊缆系的下压流阻范围为0～10 N,由于整个锚泊系统总净浮力约6 600 N,故垂向海流对于缆绳整体姿态影响程度较小。

根据图中主浮体深度曲线可之,4月19日锚泊系统下放至姿态稳定后,主浮体位置与预设水下400 m基本相同。锚泊系统布放期间有3个沉降峰值,分别为(1)2017年7月5日下潜至水下577 m,实际下沉177 m,受东北方向海流影响; (2)2017年8月28日下潜至水下635 m,实际下沉235 m,受西南方向海流影响; (3)2018年4月27日下潜至水下829 m,实际下沉429 m,受西南方向海流影响。

2.2 潜标布放位置海水密度变化影响

全布放时间段200 m至500 m SBE37CTD解析数据结果图如图7所示,由盐度和水温与深度对比曲线可以得知环境水温均低于17 ℃,代入相应海水盐度公式。

图7 200～500 m SBE37CTD温度盐度与深度对比图Fig.7 Charts comparing the 200～500 m SBE37CTD temperature salinity and depth values

其中,S为盐度;Ω为海水比重;T为水温。计算可得200 m至500 m系泊大深度沉降时海 水密度变化范围为1 023.5 kg/m3至 1 026.6 kg/m3。

受海水密度影响的因素包括: a) 根据公式(1)、(3)、(4)可知,海水密度会对流阻情况产生影响,流阻系数取1.2时,根据实际数据可知在2018年4月27日主浮体受流阻最大,U和V方向流阻总变化范围分别为66.8～67.0 N 和 883.9～886.5 N,波动区间较小; b) 根据阿基米德浮力公式知,海水密度的变化导致锚泊浮力改变,结合0～1 000 m深海水密度变化数据,根据实际布放潜标锚系结构体积计算,可得知总浮力的全时段变化范围为-50～50 N,在特定涡旋环境下变化范围为0～50 N,已知锚泊系统总净浮力约为6 600 N,在有中尺度涡的海流环境影响下,50～100 N的浮力改变足以改变本系泊的沉降深度约10 m左右。

2.3 沉降异常深日期表征海流影响

2017年7月5日、8月28日以及2018年4月27日海面高度与流速图[28]如图8、图9、图10所示。

图8 2017 年 7 月 5 日海面高度与流速图(25°～30°N,130°～150°E)Fig.8 Map of sea surface height and flow velocity on July 5,2017 (25°～30°N,130°～150°E)

图9 2017 年 8 月 28 日海面高度与流速图(25°～30°N,130°～150°E)Fig.9 Map of the sea surface height and flow velocity on August 28,2017 (25°～30°N,130°～150°E)

图10 2018 年 4 月 23 日海面高度与流速图(25°～30°N,130°～150°E)Fig.10 Map of the sea surface height and flow velocity on April 23,2018 (25°～30°N,130°～150°E)

在潜标布放海域涡旋的移动方向大部分为由东至西,偶有南北移动方向中尺度涡的形成。系泊在7月5日、8月28日以及次年4月23日,分别处于冷涡和暖涡中流速极大位置,且位于迁移方向一侧,其合成海流流速高达2～3 m/s。

为证明系泊沉降姿态确与中尺度涡迁移引起流速变化相关,拟通过对模型施加海流状态,得出其姿态模拟结果并与实际深海设备数据记录情况进行对比验证。

3 深海潜标锚泊系统仿真与验证

3.1 系统模型的建立

利用ADAMS软件建立1∶1真实锚泊模型,其水中无流情况下模型如图11所示。

图11 深海潜标锚泊系统1∶1模型图Fig.11 1∶1 Model of the deep-sea submersible anchor mooring system

其中水流阻力施加的参数包括 a) 海水密度1 025 kg/m3; b) 各个结构流阻系数施加情况如表3所示。

表3 潜标各结构雷诺数变化范围及阻力系数Tab.3 Reynolds number variation range and resistance coefficient of each structure of submersible standard

模型初始位置为前一天的夹角姿态,边界条件采取海流载荷数据,并在模型中以流阻形式直接施加于相应模型上。施加后xy二维方向和zy二维方向顶端玻璃浮球模型及受力位置如图12所示。

图12 水阻施加后顶端模型图Fig.12 Top model diagrams after the application of water resistance

3.2 锚泊系统缆系张力结果与修正

施加 3个异常沉降位的海流环境,进行系泊沉降模拟(流阻系数=1.2),2018年4月27日模型沉降最深、张力最大。比较设备与缆绳节点受力情况进行后摘取该日结果(1) 1 800 m处玻璃浮球底部节点(灰色虚线); (2) 3 000 m处玻璃浮球组底部节点(玫红色虚线); (3) 4 000 m处玻璃浮球组底部节点(浅蓝色虚线); (4) 重力锚顶部节点(深蓝色虚线),4个点位张力曲线,如图13所示。

图13 2018年4月27日锚系姿态仿真缆绳张力峰值图Fig.13 Diagram of the cable tension peaks of the anchoring attitude simulation on April 27,2018

剔除软件计算器运行方式导致的奇异值和运行程序值。锚系张力最大处为3 000 m处玻璃浮球底部节点和重力锚顶部节点,其张力大小多次达到 0.7 t左右,使用的凯夫拉缆绳破断力为 6.0 t以上,水中拉伸系数小于等于3.3%,500 kg拉力下拉伸系数小于等于0.78%。故6 000 m缆绳拉伸量最大为83.4 m,将其代入修正设备位置进行下步运算。

3.3 流阻系数值的确定

筛选两个系泊沉降较小日对比全布放时间系统沉降最深日,分别为1、2、3号模型,施加相应海流环境,分别代入5组不同流阻系数C,将主浮体深度位置模拟结果统计如表4。

由模拟结果可知,1号、2号模型对比实际沉降深度皆指向C位于1.1和1.2之间,依据文献资料[29]可查C在1.1和1.2之间。同时从表4和图14中可知,3号模型C的选取更倾向于C=1.3,且流阻系数对沉降深度和横摇的影响效果更明显。故推断,在海流流速更大的位置,若代入ADCP实测海流值应使用更大流阻系数,流阻系数的正确选择更有利于推测锚泊系统沉降的位置。

表4 不同流阻系数对浮球位置深度影响对比表Tab.4 Comparison of the influence of different flow resistance coefficients on float position depth

图14为施加不同流阻系数后,模拟缆绳的倾斜情况,其海流数据依据为2018年4月27日二维海流环境。

图14 不同流阻系数影响下浮球位置深度对比图(2018年4月27日环境载荷)Fig.14 Comparison of the position and depth of the float under the influence of different flow resistance coefficients (Environmental load on April 27,2018)

3.4 不同环境载荷的模型模拟结果与原因探究

在潜标锚系系统布放成功的1年2个月时间内,3个中尺度涡影响沉降深度时长约为98 d。为验证锚系系统在海流流正常的情况下模型计算准确度,拟选取每月第一天的流速作为参考,计算二维环境和三维环境海流参数的计算机模拟下,其 200 m、300 m、400 m、500 m位置的CTD沉降位置,并与实际沉降位置情况相对比,对比情况图15所示。

图15 三维海流环境模型中CTD沉降计算情况对比图Fig.15 Comparison of the CTD settlements calculated by 3D ocean current environmental models

根据数据计算可知,200 m、300 m、400 m和500 m的 CTD三维环境模型沉降模拟误差分别为9.27 m、8.01 m、10.33 m和7.90 m,不同位置误差相似原因可能由于CTD绑定位置紧凑。误差范围为-11～11 m,最大误差长度占模型总长0.18%。

随后,将 3个大深度沉降日期的三维海流环境施加于模型,其主浮体沉降情况比较表如表5所示,CTD模型与实际沉降深度对比如图16所示,模型偏移情况结果图如图17所示。

图16 200～500 m SBR37CTD实际沉降深度与潜标系统模型沉降情况对比图Fig.16 Comparison of the actual settlement depths of 200～500 m SBR37CTD with those obtained by the submersible target system model

图17 三个海流异常日潜标锚泊模型模拟锚泊主浮体横摇结果对比图Fig.17 Comparison of the simulated mooring results of the two-dimensional and three-dimensional ocean current environment for three outlier days

表5 400 m主浮体实际沉降深度与潜标系统模型沉降情况对比表Tab.5 Comparison of the actual settlement depths of the 400-m main floating body with those obtained by the submersible target system model

沉降深度异常日主浮体和200～500 m 4个CTD在三维海流环境模型中沉降模拟平均误差为 15 m,相对于海流较小的海流环境参数,海流异常大的情况下模拟误差明显增加,增加峰值前后两日的海流环境模拟计算,经过统计可知潜标锚泊系统模型在施加 ADCP真实海流后计算误差范围为-25～25 m,最大误差长度占模型总长0.41%。

根据 3个异常值位置的计算可发现,在深海缆系锚泊系统位于涡旋边缘尤其是位于迎涡方向时,由于海面至海深1 000 m海流流速显著增大,潜标锚系系统在当日最大流速情况下沉降至实际沉降深度,其沉降情况与前一日海流流速基本无关(忽略垂向流和海水密度变化影响)。

4 结论

综合上述分析可以得出以下结论:

1) 采用ADAMS进行6 000 m潜标锚泊系统沉降模拟,通过运用柔性体有限元方式编制模型,施加以数据为支撑的三维海流分布环境,证明在没有较强海流和中尺度涡的影响下,锚泊潜标系统模拟误差介于-11 m至11 m之间,此时流阻系数选择1.1～1.2; 当海流明显增大,海流流速(单向)在 200 m 至800 m仍能达到0.4 m/s时,考虑系统是否位于中尺度涡附近,其分析误差介于-25 m至25 m之间,为缆绳总长的 0.41%,此时流阻系数选择 1.3。此建模方式可应用于多种锚泊布放阶段系统姿态模拟,如潜标布放拖行阶段、重力锚入水后沉至海底阶段等。

2) 深水潜标锚泊系统计算误差原因可能存在于:(a) 软件模拟海流环境时忽略了 1 200 m以下海流(海流计数据受缆绳波动速度影响); (b) 缆绳实际各向异性的材料性质与模拟模型的差异; (c) 浮力调整的误差。

3) 影响系泊沉降最大因素为当日水平方向海流,当其位于漩涡边缘时,依据流速的大小,将会出现相应沉降异常值,模型对张力的计算将对缆系锚泊系统各个结构的耐压性能和连接部件的抗拉性能提出更高要求。大洋中尺度涡的迁移位置预测可有效提高深水潜标寿命,规避风险,精确预测下潜深度,更高效地获得海洋信息。