随着科学技术的发展和人类对油气资源的需求不断增加，海洋油气资源的开发逐渐向深水延伸，海洋浮式结构物的定位技术面临巨大的挑战。传统的悬链式锚泊系统通常采用三段浮容重、刚度和长度均不相同的钢链和钢缆组合而成的复合式锚泊线，利用钢链和钢缆的自身重量为上部浮体提供恢复力。但随着海洋平台工作水深的增加，浮式结构物需要更长的锚泊线系泊，从而导致更大的系泊半径、较低的回复效率、更大的缆索自重和较小的平台有效承载能力等缺陷。锚泊系统的造价更高，安装难度也越来越大，这限制了传统的悬链式锚泊系统在深水和超深水中的应用^[1]。新型人工合成纤维缆，因其自身重量更轻、成本更低，且具有较高的断裂强度，能够大大减少锚泊系统的自身重量、提升锚泊系统的力学性能，在深水和超深水应用中具有明显的优势而倍受关注。目前，合成纤维缆已被广泛应用于深水平台的锚泊系统中，用于代替钢链-钢缆-钢链复合式锚泊线中的钢缆，并被证明具有良好的经济性能^[2]。

目前，无论在理论、实验和数值分析方面，针对这种新型深水复合式锚泊系统的研究都处于起步阶段。

1 合成纤维缆动刚度特性及求解方法

复合式锚泊线中的钢缆和钢链在动力学计算中可假定为线弹性材料，而合成纤维缆却具有复杂的材料非线性特性。纤维缆材料的非线性表现为弹性模量的非定值，它随着系缆的平均张力、张力变化幅值和变化周期等因素的变化而变化，这使得准确把握纤维缆在复杂海洋环境条件下的动力响应非常困难。随着合成纤维缆在海洋工程锚泊系统中的应用越来越广泛，对纤维缆材料的研究也越来越受到重视，人们开始尝试将系缆的材料特性应用到锚泊线的动力学分析中。

Del Vacchio^[3]对聚酯纤维缆进行了模型实验，指出平均张力、张力幅值和载荷周期是影响弹性模量的主要因素，并给出了常温、循环载荷作用下纤维缆的弹性模量经验公式。Bosman等^[4]以Del Vacchio的经验公式为基础，通过实验研究发现平均载荷是影响动刚度的主要因素。Casey^[5-6]对聚酯纤维缆进行了大量的模型实验，给出了多组经验公式参数。Kim^[7]采用迭代法对聚酯纤维缆的动刚度进行求解，并将动刚度特性引入到系缆的动力响应分析中。Tahar等^[8]在柔性杆理论和有限元法的基础上，引入系缆的动刚度特性对聚酯纤维系缆进行了分析。Francois等^[9]通过模型试验指出在随机载荷作用下，平均张力是影响动刚度值的最主要因素。刘海笑和黄泽伟等^[10-11]、张火明等^[12]在对绷紧式系泊系统的数值分析中考虑了纤维系缆的动刚度特性，改进了系缆动刚度和动张力的计算方法，并以同一座工作于1 500 m水深的Spar平台为例进行了计算。

1.1 合成纤维系缆的动刚度特性

对于复合式锚泊线中的合成纤维缆部分，其响应和性能主要取决于系缆的轴向刚度(E×A) 特性(其中E为系缆的弹性模量，A为横截面积)，轴向刚度的准确表达是精确计算平台运动及系缆张力响应的基础。通常引入一个经验公式来描述合成纤维缆在循环载荷作用下的动刚度特性。

$ k=\frac{EA}{\text{MBL}}={\alpha }''+{\beta }''{{L}_{m}}-{\gamma }''{{L}_{a}} $

(1)

式中：L_m为平均张力占最小破断强度(MBL) 的百分比；L_a为动态张力幅值占MBL的百分比；系数α"、β"和γ"是与合成纤维缆构造相关的参数。

1.2 传统动刚度求解

纤维缆动刚度的求解分两步进行：

1) 求平均张力。用静刚度模型k₀=(E×A)₀/MBL (k₀是由制造商提供的定值，一般在较低的单向载荷下测得，与缆绳自身特性相关) 计算平均张力，得到L_m。若浮体在平衡位置左右做简谐振动，取系缆的初始预张力即为平均张力。

2) 迭代计算动刚度。由静刚度模型可以计算出系缆平均张力，也可求得张力变化幅值L_a，但该L_a与真实动刚度情况下的L_a相差较大，须通过迭代的方法求解，将L_m和L_a代入式(1) 中，计算得到一个新的动刚度值k₁。利用k₁计算得到的新的系缆张力及L_a，再将L_a代入式(1)，求解动刚度值k₂…，如此反复迭代，直到相邻两次的动刚度计算结果k_i-k_i-1小于设定的容差值，认为计算收敛，k_i即为所求得的动刚度值。

1.3 改进的动刚度求解方法

传统动刚度计算^[10-12]中将整根纤维系缆作为一条直线考虑，假定纤维缆上的应力和应变处处相同；认为合成纤维缆的密度与海水密度十分接近，因而不考虑缆的自重，且系缆所受的流体动力相对于轴向张力较小，亦忽略不计；计算中通常假设上部浮体在平衡位置附近做简谐振动，将系缆的初始预张力作为平均张力。传统动刚度求解方法虽然大大简化了计算，但也损失了计算精度。

首先，尽管合成纤维缆自身重量较轻，但深水系泊中的缆长通常较长，缆索自重产生的总体作用仍然较大。以表 1中的聚酯纤维缆为例，其总浮重可达0.085 26×2 000.0=170.52 kN，为预张力的(2 308 kN) 的7.39%，忽略这部分的影响显然是不合适的；其次，纤维缆在外界激励下的形状并不是一条直线，系缆上的张力同一时刻也并非处处相等，下部点的响应比其上部点要滞后；最后，传统方法无法考察流体动力对缆索的动力响应。

另外，传统动刚度计算方法仅适用于由合成纤维缆组成的单成份锚泊线的计算，对于形式为锚链-纤维缆-锚链的复合式锚泊线，传统动刚度计算方法不再适用。为精确计算复合式锚泊线上的动张力，应将合成纤维缆的动刚度特性引入到锚泊线的动力学分析中。

本文采用分段动刚度的方法，将合成纤维缆按与钢链相同的方法进行空间离散，建立纤维缆段的动力学模型，对每一缆段采用误差控制的迭代方法求解其动刚度和动张力；在动力学模型中充分考虑流体动力、自身重量和海流等因素的作用；基于统计的方法计算平均张力，即记录每一缆段在过去一段时间的动张力，取该记录的平均值作为平均张力，并将其用于下一时刻该缆段动刚度的计算；为减少迭代次数，以缆段上一时刻的动刚度值作为下一时刻迭代的初始值。

2 锚泊线动力学模型

锚泊系统的动力学建模方法有集中质量法、有限元法和有限差分法等。其中，集中质量法^[13-14]因物理意义明确，算法简单易懂，具有广泛的适用性及扩展性而得到广泛应用。Chai等^[15]将集中质量法进行了扩展，并将弯矩、扭矩、与海底的接触问题等加入到海洋缆索的计算模型中。王飞^[16-18]和朱克强等^[19-20]基于集中质量法建立了海洋缆索的动力学模型，同时考虑了弯矩、拖缆-海底接触等的响应，实现了缆索收放过程的模拟。

2.1 坐标系

锚泊系统的动力学模型应该建立在合适的坐标系下，根据需要建立如图 1所示的惯性坐标系o-xyz和局部坐标系i-btn。惯性坐标系是空间固定的坐标系，所有的计算均转换到该坐标系下进行，其原点位于锚泊系统末端与锚的连接处，长度记为s=0。局部坐标系附在锚泊线上，t轴为微元的切线方向，指向长度s增加方向，n和b分别为法向和副法向；欧拉角(φ, θ) 为微元段的姿态角。两个坐标系均为右手系，通过姿态角(φ, θ) 进行关联。局部坐标系到惯性坐标系的转换矩阵为：

$ \mathit{\boldsymbol{A=}}\left[\begin{matrix} \cos \ \theta &\sin \ \theta \cos \ \varphi &-\sin \ \theta sin\ \varphi \\ -\sin \ \theta &\cos \ \theta \cos \ \varphi &-\cos \ \theta sin\ \varphi \\ 0&sin\ \varphi &\cos \ \varphi \\ \end{matrix} \right] $

式中：[x y z]^T=A[b t n]^T，矩阵A为单位正交矩阵，其逆矩阵为其转置矩阵。

2.2 动力学模型

为建立锚泊系统的动力学模型，将锚泊线在空间上离散为一系列节点。锚泊线总长度为S，末端s=0为第i=0个节点，上端点处s=S，为第i=N个节点。任取一微元段ds进行受力分析并应用牛顿第二定律，得到第i个节点的控制方程：

$ {{M}_{i}}{{{\ddot{x}}}_{i}}=\Delta {{T}_{i}}+{{B}_{i}}+{{G}_{i}}+{{F}_{D\ i}}+{{F}_{bfi}} $

(2)

式中：${{{\ddot{x}}}_{i}}$为第i个节点位置对时间的二阶导数，即加速度；ΔT_i为节点受到的张力；B_i和G_i分别为第i个节点的浮力和重力；F_Di为节点受到的流体动力；F_bfi为海底的作用力；M_i为节点的质量矩阵。

2.2.1 质量矩阵

$ {{M}_{i}}={{m}_{i}}I+{{M}_{ai}} $

(3)

式中：惯性质量m_i=(μ_i-1/2l_i-1/2+μ_i+1/2l_i+1/2)/2；附加质量：$\begin{align} &{{M}_{ai}}=\left( {{M}_{ai-1/2}}+{{M}_{ai+1/2}} \right)/2= \\ &\left( \rho {{k}_{i-1/2}}{{l}_{i-1/2}}\cdot {{P}_{i-1/2}}+\rho {{k}_{i+1/2}}{{l}_{i+1/2}}{{\sigma }_{i+1/2}}\cdot {{P}_{i+1/2}} \right)/2; \\ &\ \ \ \ \mathit{\boldsymbol{P=}} \\ &\left[\begin{matrix} 1-{{\sin }^{2}}\theta {{\cos }^{2}}\varphi &-\sin \ \theta \cos \ \theta {{\cos }^{2}}\varphi &-\sin \ \theta \sin \ \varphi \cos \ \varphi \\ -\sin \ \theta \cos \ \theta {{\cos }^{2}}\varphi &1-{{\cos }^{2}}\theta {{\cos }^{2}}\varphi &-\cos \ \theta \sin \ \varphi \cos \ \varphi \\ -\sin \ \theta \sin \ \varphi \cos \ \varphi &-\cos \ \theta \sin \ \varphi \cos \ \varphi &{{\cos }^{2}}\varphi \\ \end{matrix} \right] \\ \end{align}$；μ，l，ρ，k，σ分别为系缆单位长度质量，节点间的长度，流体密度，附加质量系数和横截面积；下标i-1/2和i+1/2分别表示节点i-1到i、节点i到i-1之间的物理量。

2.2.2 浮力和重力

$ {{B}_{i}}+{{G}_{i}}=\rho g\left( {{l}_{i-1/2}}{{\sigma }_{i-1/2}}+{{l}_{i+1/2}}{{\sigma }_{i+1/2}} \right)/2-{{m}_{i}}g $

(4)

2.2.3 流体动力

按Ablow^[21]和Huang^[22]的方法：

$ {{F}_{D\ i}}=F_{D}^{\ \ \ i-1/2}+F_{D}^{\ \ \ i+1/2} $

(5)

在局部坐标系下，节点的流体动力为$\left\{ \begin{align} &{{F}_{bi+1/2}}=-\frac{1}{2}\rho {{C}_{b}}{{d}_{i+1/2}}{{l}_{i+1/2}}\left| {{\left( V-U \right)}_{b}} \right|{{\left( V-U \right)}_{b}} \\ &{{F}_{ti+1/2}}=-\frac{1}{2}\rho \text{ }\!\!\pi\!\!\text{ }{{C}_{t}}{{d}_{i+1/2}}{{l}_{i+1/2}}\left| {{\left( V-U \right)}_{t}} \right|{{\left( V-U \right)}_{t}} \\ &{{F}_{ni+1/2}}=-\frac{1}{2}\rho {{C}_{n}}{{d}_{i+1/2}}{{l}_{i+1/2}}\left| {{\left( V-U \right)}_{n}} \right|{{\left( V-U \right)}_{n}} \\ \end{align} \right.$式中：[(V-U)_b (V-U)_t (V-U)_n]^T=A^-1(V-U)，(V-U) 为在惯性坐标系下节点相对于周围流体的速度，V为节点的运动速度，U为流体速度。将局部坐标系下的流体动力转换到惯性坐标系下：

$ \mathit{\boldsymbol{F}}_{D}^{\ i+1/2}=\mathit{\boldsymbol{A}}{{\left[{{F}_{bi+1/2}}\ \ {{F}_{ti+1/2}}\ \ {{F}_{ni+1/2}} \right]}^{\rm{T}}} $

2.2.4 张力

张力由链或缆的材料特性和形变决定。合成纤维系缆应力-应变关系不是简单的线性关系，它受系缆构造形式、材料、载荷类型等影响。为将纤维缆的材料特性引入到锚泊线的动力分析中，采用误差控制的迭代算法求解微元的动刚度和动张力：

$ \left\{ \begin{align} &\Delta T={{T}_{i+1/2}}-{{T}_{i-1/2}} \\ &{{\left( E\sigma \right)}_{i+1/2}}={{k}_{i+1/2}}\times \text{MBL, }{{\left( E\sigma \right)}_{i+1/2}}={{k}_{i-1/2}}\times \text{MBL} \\ &{{T}_{i+1/2}}={{\left( E\sigma \right)}_{i+1/2}}\times {{\varepsilon }_{i+1/2}}, {{T}_{i-1/2}}={{\left( E\sigma \right)}_{i-1/2}}\times {{\varepsilon }_{i-1/2}} \\ &{{\varepsilon }_{i+1/2}}= \\ &^{\underline{\sqrt{{{\left( {{x}_{i+1}}-{{x}_{i}} \right)}^{2}}+{{\left( {{y}_{i+1}}-{{y}_{i}} \right)}^{2}}+{{\left( {{z}_{i+1}}-{{z}_{i}} \right)}^{2}}}}}/{{l}_{i+1/2}}-1 \\ \end{align} \right. $

(6)

钢缆和钢链假定为线弹性材料，其本构关系采用虎克定律：

$ \left\{ \begin{align} &\Delta T={{T}_{i+1/2}}-{{T}_{i-1/2}} \\ &{{T}_{i+1/2}}={{\sigma }_{i+1/2}}E{{\varepsilon }_{i+1/2}}, {{T}_{i-1/2}}={{\sigma }_{i-1/2}}E{{\varepsilon }_{i-1/2}} \\ &{{\varepsilon }_{i+1/2}}= \\ &\sqrt{{{\left( {{x}_{i+1}}-{{x}_{i}} \right)}^{2}}+{{\left( {{y}_{i+1}}-{{y}_{i}} \right)}^{2}}+{{\left( {{z}_{i+1}}-{{z}_{i}} \right)}^{2}}}/{{l}_{i+1/2}}-1 \\ \end{align} \right. $

(7)

2.2.5 与海底的相互作用

采用常洪波的方法^[23]，将海底土壤视为线性的弹性基础，用库仑定律描述摩擦力。节点i的垂向坐标为z_i，该处海底坐标为z_bi，当z_i < z_bi时，海底对节点的支持力F_bi为：

$ {{F}_{bi}}=k\left( {{z}_{bi}}-{{z}_{i}} \right)\cdot n+c{{{\dot{x}}}_{i}}n $

(8)

摩擦力的计算公式为：

$ {{F}_{fi}}=\left\{ \begin{align} &-\mu {{F}_{bi}}\frac{{{v}_{ti}}}{\left| {{v}_{ti}} \right|}{{v}_{ti}}>{{v}_{\text{lim}}} \\ &-\mu {{F}_{bi}}\frac{{{v}_{ti}}}{\left| {{v}_{\text{lim}}} \right|}{{v}_{ti}}\le {{v}_{\text{lim}}} \\ \end{align} \right. $

(9)

式中：k为海底的等效刚度系数，c为等效阻尼系数，n为海底的外法线方向，μ为锚泊线与海床面的摩擦系数，v_ti为与海底接触的节点的切向速度，设置临界速度v_lim应尽可能小。

2.3 动力学模型数值解算

1) 边界条件。锚端边界条件设置为固定端；锚泊线顶端的位置和速度与平台的运动保持一致：

$ \left\{ \begin{array}{l} {x_N} = {x_s}\left( t \right)\\ {y_N} = {y_s}\left( t \right)\\ {z_N} = {z_s}\left( t \right) \end{array} \right.\;\;\;\;\;{\rm{and}}\;\;\;\;\;\left\{ \begin{array}{l} {{\dot x}_N} = {u_s}\left( t \right)\\ {{\dot y}_N} = {v_s}\left( t \right)\\ {{\dot z}_N} = {w_s}\left( t \right) \end{array} \right. $

(10)

式中：x_s、y_s、z_s、u_s、v_s、w_s分别为锚泊线顶端的位置和速度，它们是时间的函数。

2) 初始条件。确定节点初始时刻位置和速度：

$ \left\{ \begin{array}{l} {x_i}\left( 0 \right) = x_i^0\\ {y_i}\left( 0 \right) = y_i^0\\ {z_i}\left( 0 \right) = z_i^0 \end{array} \right.\;\;\;\;\;{\rm{and}}\;\;\;\;\;\left\{ \begin{array}{l} {{\dot x}_i}\left( 0 \right) = \dot x_i^0\\ {{\dot y}_i}\left( 0 \right) = \dot y_i^0\\ {{\dot z}_i}\left( 0 \right) = \dot z_i^0 \end{array} \right. $

(11)

式中：等式右侧部分为给定的初始值。

3) 数值求解。联立控制方程(2)，再加上初始条件，并由v=dx/dt，得到完整的偏微分方程组：

$ \left\{ \begin{align} &\frac{\text{d}{{{\dot{x}}}_{i}}}{\text{d}t}=M_{ai}^{-1}\cdot {{F}_{i}} \\ &\frac{\text{d}{{x}_{i}}}{\text{d}t}={{{\dot{x}}}_{i}} \\ \end{align} \right. $

(12)

对方程组采用四阶龙格-库塔法求解，由各节点在t_n时刻运动状态即可得到t_n+1=t_n+Δt时刻的运动状态。

3 计算案例

以一座工作水深为1 500 m的海洋平台为例进行计算，该平台锚泊系统采用复合式锚泊线。锚泊线参数如表 1所示，顶端的初始预张力为2 308 kN，锚泊线结构如图 2所示。聚酯纤维缆的材料特性参数α"=14.469，β"=0.211 3，γ"=0.269 7，准静刚度值k₀=12.2^[2]。

将锚泊线按每段25.0 m离散为158段，时间步长取2 ms。将k₀=12.2代入到动力学模型中的纤维缆段进行计算，复合式锚泊线在初始预张力作用下的构形如图 2所示。由于聚酯纤维缆段重量较轻，其空间形状接近为一条直线。将锚泊线首端位置固定，可得到静刚度下锚泊线首尾两端的张力-时间历程如图 3所示。

若采用改进的动刚度方法对锚泊线进行动力学计算，锚泊线首尾两端的张力-时间历程如图 4和图 5所示，锚泊线首尾两端的动张力比静刚度下计算结果分别大5.50%和7.12%；纤维缆段上下两端的动张力计算结果如图 6和图 7所示，比静刚度下的计算结果分别大6.76%和7.13%。

图 8为基于统计的方法得到的纤维缆段上的平均张力，该值沿缆长的方向增加。图 9为平均张力占最小破断强度(MBL) 的百分比(L_m)，该值介于13.47~14.06。在合成纤维缆的动力特性分析中，当L_m＞10.0时就必须考虑纤维缆的动刚度特性。

图 10为纤维缆两端的张力变化幅值占最小破断强度(MBL) 的百分比(L_a)，该值在0值附近振动。L_m和L_a决定了不同缆段上的动刚度计算结果。图 11为纤维缆段两端的动刚度值，分别位于17.44和17.315附近，均比静刚度值k₀=12.2大得多。

假设平台沿x方向运动，锚泊线上端随平台主体发生位移，位移随着时间的变化历程为正弦函数为x(t)=x₀sin (2πt/T)。x₀取5.0 m，T取10 s。

图 12和图 13分别为施加激励后锚泊线两端动张力的计算结果及其与静刚度下计算结果的比较情况。锚泊线首端最大和最小张力分别比静刚度下大5.79%和7.91%，尾端最大和最小张力分别比静刚度下大12.34%和7.26%。

图 14和图 15分别为纤维缆两端动张力计算结果，纤维缆上端最大和最小张力分别比静刚度下大12.32%和4.74%，下端最大和最小张力分别比静刚度下大7.82%和6.70%。

图 16和图 17分别为纤维缆两端的平均张力和张力幅值占最小破断强度(MBL) 的百分比，二者构成影响动刚度计算结果的主要因素。

图 18为纤维缆两端的动刚度值的计算结果，可以看出动刚度值比静刚度值大得多，且随着上端激励的变化而周期变化。图 14、15、17和18表明纤维缆段下端的动张力变化幅值和动刚度变化幅值均比上端大。

4 结论

本文建立了复合式锚泊系统的动力学模型，将合成纤维缆的动刚度特性应用到组合式锚泊线的数值计算中，通过对一座工作水深为1 500 m的组合式锚泊线进行计算，得到以下结论：

1) 传统动刚度计算中因为存在大量的假设和忽略，降低了计算精度，且不适用于不均匀缆和复合式锚泊线中纤维缆的动力学计算。为了精确计算复合式锚泊线上的动张力，应该建立复合锚泊线的动力学模型，并将合成纤维缆的动刚度特性引入到锚泊线的动力学分析中。

2) 在纤维缆的动力学计算中应采用分段动刚度的方法计算各缆段上的动刚度和动张力。传统方法中取恒定预张力作为平均张力的方法不再适用，可采用基于统计的方法计算各缆段的平均张力，并用于该缆段上动刚度的计算。

3) 采用动刚度方法得到的复合式锚泊线动张力远大于静刚度下的计算结果，因而在合成纤维缆的动力学计算中必须考虑其动刚度特性。

4) 最终计算结果显示，锚泊线首端无激励时首尾两端的动张力比采用静刚度方法时分别大5.5%和小7.12%；在正弦激励作用下锚泊线首端最大和最小张力比静刚度方法分别大5.79%和7.91%，尾端最大和最小张力比静刚度方法分别大12.34%和7.26%。

本文计算合成纤维缆动刚度和动张力的方法，更有合理性，并考虑到了各因素的影响，能适用于不均匀缆和复合式锚泊线的计算，对于新型复合式锚泊系统的数值分析和工程应用具有重要的意义。