水下纜索運(yùn)動(dòng)建模與仿真分析方法

張 斌，劉延平

（1. 中國(guó)船舶集團(tuán)有限公司第七一〇研究所，湖北 宜昌 443003；2. 清江創(chuàng)新中心，湖北 武漢 430076）

0 引言

隨著對(duì)海洋開(kāi)發(fā)利用的不斷深入，各類(lèi)水下無(wú)人錨泊駐留系統(tǒng)、海洋拖曳系統(tǒng)逐漸成為相關(guān)熱點(diǎn)研究領(lǐng)域，準(zhǔn)確預(yù)報(bào)分析錨泊、拖曳等系統(tǒng)在海流等外界環(huán)境擾動(dòng)下的動(dòng)態(tài)運(yùn)動(dòng)響應(yīng)成為其工程設(shè)計(jì)的關(guān)鍵。其中，建立用于連接設(shè)備、傳遞或緩沖載荷、傳輸能源和信號(hào)等用途的柔性纜索的運(yùn)動(dòng)學(xué)與動(dòng)力學(xué)精確建模是一項(xiàng)難點(diǎn)技術(shù)。

水下柔性纜索的力學(xué)、運(yùn)動(dòng)學(xué)特征與錨泊、拖曳系統(tǒng)整體狀態(tài)耦合相關(guān)，并交互影響，且隨著外部載荷變化呈現(xiàn)非線(xiàn)性時(shí)變狀態(tài)，難以對(duì)描述其狀態(tài)的微分方程進(jìn)行直接求解，通常在滿(mǎn)足動(dòng)力學(xué)條件和連續(xù)性條件下，采用有限差分[1-2]、集中質(zhì)量[3-4]、有限段[5-6]等空間域內(nèi)的物理離散方法對(duì)連續(xù)纜索控制方程進(jìn)行近似求解。

在水下纜索中內(nèi)部張力占主導(dǎo)作用情況下，可忽略纜索單元空間旋轉(zhuǎn)運(yùn)動(dòng)，上述線(xiàn)性離散化求解方法可取得較好的計(jì)算精度和準(zhǔn)確性。但在處理多點(diǎn)懸鏈系留、錨泊狀態(tài)振蕩急劇變化、拖船大范圍轉(zhuǎn)向等易于出現(xiàn)纜索低張力松弛、非線(xiàn)性大形變等問(wèn)題時(shí)，線(xiàn)性離散單元無(wú)法準(zhǔn)確反映描述纜索空間曲率引起的彎矩、扭矩作用，難以描述纜索真實(shí)空間狀態(tài)與內(nèi)力狀態(tài)，進(jìn)而出現(xiàn)數(shù)值求解奇異。通過(guò)選取大量離散節(jié)點(diǎn)進(jìn)行更加細(xì)的線(xiàn)性單元?jiǎng)澐謥?lái)逐漸逼近纜索真實(shí)輪廓，不僅算力消耗劇增，也面臨數(shù)值積分截?cái)嗾`差與舍入誤差累加問(wèn)題導(dǎo)致的計(jì)算錯(cuò)誤。

因此，本文改進(jìn)柔性纜索的線(xiàn)性離散化單元為高階非線(xiàn)性單元，采用三次樣條插值描述纜索微元的空間形態(tài)，通過(guò)加權(quán)余量法解決插值方法僅能在離散節(jié)點(diǎn)位置滿(mǎn)足纜索狀態(tài)微分方程問(wèn)題，消除微元長(zhǎng)度范圍內(nèi)殘差，提供一種高效、準(zhǔn)確求解水下纜索復(fù)雜運(yùn)動(dòng)狀態(tài)的新方法。

1 水下纜索空間運(yùn)動(dòng)數(shù)學(xué)模型

引入圖1中所示的慣性坐標(biāo)系O（xe，ye，ze）、纜索局部坐標(biāo)系 S（q1，q2，q3）以及弗萊納坐標(biāo)系 S（t，n，b）。

圖1 坐標(biāo)系選取及水下纜索微元力學(xué)狀態(tài)Fig.1 Coordinate systems involved and micro-elements force of the cable element

對(duì)于長(zhǎng)度為ds的纜索微元，其在流體作用力、自身體積力、內(nèi)部張力等共同作用下，大地坐標(biāo)系中空間平移與旋轉(zhuǎn)運(yùn)動(dòng)的微分控制方程如下[7]：

式中：Fc為微元內(nèi)部張力矢量；h為微元長(zhǎng)度上的流體作用力；wc為微元水下濕重矢量；f為其他外部作用力；MI為纜索微元質(zhì)量矩陣；r（s，t）為纜索微元的空間構(gòu)型曲線(xiàn)矢量；Nc為微元內(nèi)部的彎矩、扭矩矢量；ω為纜索微元的旋轉(zhuǎn)角速度矢量陣。

應(yīng)綜合考慮坐標(biāo)系S（q1，q2，q3）以及S（t，n，b）本身不重合導(dǎo)致的角間距 δ（s，t）偏差，與纜索微元形狀扭轉(zhuǎn)變形共同組成總撓率[8]：

纜索微元內(nèi)部任一點(diǎn)處力矩矢量可表述為沿弗萊納坐標(biāo)系 b軸方向的彎矩以及沿 t軸方向的扭矩[9]：

式中：I為纜索橫截面慣性矩；Ip為纜索橫截面極慣性矩；G為纜索的材料剪切模量。

水下纜索大長(zhǎng)細(xì)特性使得其轉(zhuǎn)動(dòng)慣量矩陣 J本身較小，且其法向阻力系數(shù)一般以量級(jí)差距遠(yuǎn)大于切向阻力系數(shù)，進(jìn)一步限制其產(chǎn)生旋轉(zhuǎn)角加速度，大地坐標(biāo)系中其旋轉(zhuǎn)控制方程等式右半部分為2個(gè)極小值乘積，可認(rèn)為水下纜索微元內(nèi)部力與力矩形成近似平衡狀態(tài)，可得：

與弗萊納坐標(biāo)系單位矢量t叉乘，去除上述旋轉(zhuǎn)運(yùn)動(dòng)控制方程中沿切向t并對(duì)內(nèi)部力矩?zé)o影響分量部分，可得：

通過(guò)叉乘反交換率與拉格朗日公式計(jì)算上式中連續(xù)叉乘，可得到沿水下纜索微元軸向任意位置處的內(nèi)部張力矢量：

式中，TC為纜索內(nèi)部軸向拉力。

通過(guò)水下纜索微元空間構(gòu)型曲線(xiàn)r(s, t)描述弗萊納坐標(biāo)系各軸線(xiàn)矢量方向，具體形式如下：

可將內(nèi)部力 FC從弗萊納坐標(biāo)系中轉(zhuǎn)換至大地坐標(biāo)系中，簡(jiǎn)化對(duì)于彎曲、扭轉(zhuǎn)等力矩作用效果的求解。

水下纜索微元大地坐標(biāo)系中空間平移運(yùn)動(dòng)方程可寫(xiě)為如下形式，包含了微元中軸向拉力、扭轉(zhuǎn)應(yīng)力、彎曲應(yīng)力對(duì)于纜索空間形態(tài)r(s, t)的影響：

由于水下纜索中扭矩沿軸向連續(xù)且恒定，可知旋轉(zhuǎn)運(yùn)動(dòng)控制方程在弗萊納坐標(biāo)系t軸方向投影必定為0：

上式給出了坐標(biāo)系S（q1，q2，q3）以及S（t，n，b）本身不重合導(dǎo)致的角間距 δ（s，t）偏差的計(jì)算方法，進(jìn)而可以確定坐標(biāo)系 S（q1，q2，q3）的坐標(biāo)軸指向。

2 水下纜索運(yùn)動(dòng)控制方程離散化求解

2.1 水下纜索非線(xiàn)性離散化單元構(gòu)造

由連續(xù)纜索運(yùn)動(dòng)控制方程計(jì)算推導(dǎo)過(guò)程可知，坐標(biāo)系矢量軸方向、纜索形變量均可通過(guò)纜索空間形態(tài)r（s，t）及角間距δ（s，t）求解得到。

在空間域內(nèi)對(duì)水下纜索進(jìn)行離散化處理，將其劃分為N段長(zhǎng)度為ds的非線(xiàn)性高階單元。根據(jù)[A1]大地坐標(biāo)系下節(jié)點(diǎn)坐標(biāo)值 r（i）（i=1，2…N，N+1）通過(guò)3次樣條差值方法構(gòu)造試函數(shù)逼近空間樣條曲線(xiàn)的真實(shí)分布狀態(tài)，二階連續(xù)可導(dǎo)既滿(mǎn)足空間曲率連續(xù)要求，又能避免更高階函數(shù)帶來(lái)的額外方程求解問(wèn)題。

圖2 水下纜索非線(xiàn)性單元?jiǎng)澐諪ig.2 Non-linear element of the underwater cable

2.2 加權(quán)余量法求解

式中：Wvi是求解域內(nèi)的權(quán)函數(shù)；Rv為試函數(shù)偏差產(chǎn)生的計(jì)算殘數(shù)；V為去除邊界后的求解計(jì)算域。

選用Galerkin方法為求解具體方法，其權(quán)函數(shù)與基函數(shù)一致，可將微分方程轉(zhuǎn)換為一系列對(duì)稱(chēng)的代數(shù)方程，實(shí)現(xiàn)計(jì)算收斂的同時(shí)進(jìn)一步減小算力需求。對(duì)于第i（i=1，2…N，N+1）段非線(xiàn)性纜索微元，其計(jì)算域內(nèi)的余量消除方程為

受篇幅所限，纜索單元求解域內(nèi)部張力、流體力、體積力、彎曲力、扭轉(zhuǎn)力等矩陣的具體形式不再贅述[11-12]。

方程（17）-（18）完整描述了每一段纜索單元空間平移與旋轉(zhuǎn)運(yùn)動(dòng)。對(duì)于纜索節(jié)點(diǎn)i，其廣義運(yùn)動(dòng)速度和加速度在其前后兩段纜索單元運(yùn)動(dòng)控制方程中體現(xiàn)，將第i-1段、第i段單元中關(guān)于節(jié)點(diǎn)i的相關(guān)項(xiàng)疊加整合，可得節(jié)點(diǎn)1至節(jié)點(diǎn)N處運(yùn)動(dòng)方程：

式中：[H]、[W]、[B]為各離散節(jié)點(diǎn)疊加整合后的流體作用力、纜索質(zhì)量及浮力矩陣；[M]為（3N+3）階質(zhì)量矩陣；[KAXIA]、[KBEND]、[KTORQ]為各離散節(jié)點(diǎn)疊加整合后的軸向張力、彎曲作用力、扭轉(zhuǎn)作用力矩陣，為（3N+3）×（6N+6）階形式矩陣。

對(duì)于常見(jiàn)的水下錨泊系統(tǒng)、拖曳系統(tǒng)，纜索上一般無(wú)外部額外施加力矩，因此，可將非線(xiàn)性纜索單元的扭轉(zhuǎn)運(yùn)動(dòng)約束方程整合為如下形式，計(jì)算中需指定水下纜索首節(jié)點(diǎn)扭轉(zhuǎn)形變角δ（1）或末節(jié)點(diǎn)處扭轉(zhuǎn)形變角 δ（N+1）以固定一個(gè)扭轉(zhuǎn)自由度，使得該約束方程左側(cè)系數(shù)矩陣變?yōu)闈M(mǎn)秩矩陣，進(jìn)而獲得唯一解。

3 模型驗(yàn)證

3.1 Hopland拖曳試驗(yàn)驗(yàn)證

采用1992年Hopland進(jìn)行的通信光纜拖曳布放實(shí)海實(shí)驗(yàn)測(cè)試數(shù)據(jù)作為對(duì)比算例[13]，通過(guò)水下纜索運(yùn)動(dòng)模型對(duì)其重型鎧裝拖纜在拖船加減速過(guò)程中動(dòng)態(tài)變化情況。

該纜索主要物理特性如下：直徑33.2 mm，總長(zhǎng) 300 m，密度 3 121 kg/m3，彈性模量77.5×109Pa，抗彎剛度1 000 N·m2，法向阻力系數(shù)1.649，切向阻力系數(shù)0.12。

加速實(shí)驗(yàn)過(guò)程中，拖船以約1.1 kn勻航速拖動(dòng)纜索達(dá)到穩(wěn)定狀態(tài)，隨后進(jìn)入加速拖曳階段，拖船速度在60 s時(shí)間內(nèi)均勻增大至2.4 kn。計(jì)算模型仿真計(jì)算結(jié)果如圖3所示，圖中，t=0時(shí)刻的纜形圖對(duì)應(yīng)于1.1 kn均勻航速拖曳下纜索狀態(tài)，t=60 s時(shí)刻，拖船線(xiàn)性加速達(dá)到2.4 kn，由仿真結(jié)果可知，拖纜其重新趨于穩(wěn)態(tài)時(shí)間滯后約360 s。

圖3 拖船加速過(guò)程HA纜姿態(tài)變化Fig.3 Attitude change of HA cable during the towing ship accelerating process

減速實(shí)驗(yàn)過(guò)程中，拖船以約2.5 kn勻航速拖動(dòng)纜索達(dá)到穩(wěn)定狀態(tài)，隨后進(jìn)入減速拖曳階段，拖船速度在60 s時(shí)間內(nèi)均勻降低至1.0 kn。計(jì)算模型仿真計(jì)算結(jié)果如圖4所示，圖中t=0時(shí)刻的纜形圖對(duì)應(yīng)于2.5 kn均勻航速拖曳下纜索狀態(tài)，t=60 s時(shí)刻，拖船線(xiàn)性減速達(dá)到1.0 kn，由仿真結(jié)果可知，拖纜重新趨于穩(wěn)態(tài)時(shí)間滯后約460 s。

圖4 拖船減速過(guò)程HA纜姿態(tài)變化Fig.4 Attitude change of HA cable during the towing ship decelerating process

綜合上述實(shí)驗(yàn)測(cè)試?yán)|索形態(tài)變化及 Hopland實(shí)測(cè)的纜索節(jié)點(diǎn)空間位置坐標(biāo)，可以看到數(shù)值模擬數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)吻合度及變化規(guī)律高度一致，僅存在細(xì)微偏差，推測(cè)是由于實(shí)驗(yàn)測(cè)量誤差及實(shí)驗(yàn)過(guò)程中海流流速影響等導(dǎo)致。

3.2 細(xì)長(zhǎng)柔性桿件彎曲形變計(jì)算驗(yàn)證

通過(guò)與可精確求得解析值的經(jīng)典細(xì)長(zhǎng)柔性桿件受力彎曲形變模型進(jìn)行比對(duì)分析，驗(yàn)證本計(jì)算模型中高階非線(xiàn)性微元?jiǎng)澐旨凹訖?quán)余量方法消除殘差趨近物理模型真實(shí)空間形態(tài)方法的準(zhǔn)確性，也直接對(duì)模型中彎曲應(yīng)力的影響結(jié)果進(jìn)行比對(duì)分析。

本算例的物理模型為一端固定約束、一段自由的細(xì)長(zhǎng)柔性桿件，對(duì)其末端施加豎直向下集中載荷P。由材料力學(xué)相關(guān)解析求解方法可知，細(xì)長(zhǎng)桿上任一點(diǎn)處曲率半徑如下：

圖5 細(xì)長(zhǎng)柔性桿件彎曲形變計(jì)算驗(yàn)證模型Fig.5 Bending deformation validation model of flexible rod

式中：MBEND為該點(diǎn)位置彎矩；E為材料彈性模量；I為桿件的橫截面慣性矩；θ為該點(diǎn)彎曲形變后繞中心軸轉(zhuǎn)動(dòng)角度；s為撓度曲線(xiàn)弧長(zhǎng)。

上述末端節(jié)點(diǎn)的位移量、偏轉(zhuǎn)角可通過(guò)數(shù)值積分求解方法獲得準(zhǔn)確值。

本算例模型相關(guān)計(jì)算參數(shù)如下：直徑20 mm、長(zhǎng)度10 m、密度1 000 kg/m3、彈性模量2×1011Pa。采用高階非線(xiàn)性離散化模型建立桿件的運(yùn)動(dòng)控制方程，仿真分析不同集中力載荷P施加狀態(tài)下桿件形態(tài)及末端節(jié)點(diǎn)坐標(biāo)值，具體結(jié)果如圖6所示。

圖6 不同集中載荷作用下細(xì)長(zhǎng)桿平衡姿態(tài)Fig.6 Torque equilibrium attitude of the flexible rod with different loading

圖7 末端節(jié)點(diǎn)坐標(biāo)值隨時(shí)間變化情況Fig.7 Change of the end node’s coordinate value with time

圖8 桿件形態(tài)隨時(shí)間變化情況（Δt=2 s）Fig.8 Change of the attitude of the flexible rod with time（Δt=2 s）

由計(jì)算結(jié)果可知，桿件末端節(jié)點(diǎn)在集中載荷P的作用下逐漸偏轉(zhuǎn)，在彎曲應(yīng)力與載荷P共同作用下，節(jié)點(diǎn)在平衡位置附近做欠阻尼振蕩并逐漸趨于穩(wěn)定狀態(tài)。

不同計(jì)算約束條件下，使用本文建立的運(yùn)動(dòng)控制模型求解得到的末端節(jié)點(diǎn)坐標(biāo)值與解析值結(jié)果相對(duì)誤差ε如下表所示，具體計(jì)算方法為

由計(jì)算結(jié)果可知，通過(guò)高階非線(xiàn)性微元物理離散方法建立的運(yùn)動(dòng)控制方程具有較高計(jì)算精度。

表1 細(xì)長(zhǎng)柔性桿件彎曲形變計(jì)算結(jié)果比對(duì)Table 1 Comparison of the simulation value and theoretical value of the flexible rod

4 結(jié)束語(yǔ)

本文將水下纜索線(xiàn)性化離散改進(jìn)為高階非線(xiàn)性單元，建立了包含彎矩、扭矩作用的水下纜索狀態(tài)微分方程數(shù)值求解模型，通過(guò)與Hopland水下纜索拖曳試驗(yàn)數(shù)據(jù)、細(xì)長(zhǎng)柔性桿件彎曲形變計(jì)算等對(duì)比分析，驗(yàn)證了本文高階離散求解模型的準(zhǔn)確性，可為水下錨泊系留、水下拖曳等系統(tǒng)復(fù)雜運(yùn)動(dòng)狀態(tài)高效求解與工程設(shè)計(jì)提供有價(jià)值的參考。