基于模型预测控制的深海钻井立管再入井仿真分析

图1 深海钻井立管硬悬挂模式示意图

Fig.1 A hard hang-off deep-water marine riser

修正的哈密顿原理表达式为

$\int_{t_{1}}^{t_{2}}$ (δK-δV+δW)dt=0

(1)

式中:K和V分别为系统动能和势能;W为非保守力所做的虚功. K包括钻井船的动能K₁与悬挂立管(含LMRP)的动能K₂:

K₁=

$\frac{1}{2} {\overset{\cdot}{M r}}^{2}$

(2)

K₂=

$\frac{1}{2} \int_{0}^{l}$ [m+δ(z-l)M_B](

$\overset{\cdot}{w}$ (z, t)+

$\overset{\cdot}{r}$ )²dz

(3)

式中:M为钻井船质量;m为立管单位长度质量;l为立管长度.由于深海立管长径比很大,柔性很强,可采用缆绳模型进行计算^{[3⇓⇓⇓⇓⇓⇓⇓⇓-12]},所以没有弯曲应变能.立管系统势能V包括重力势能V_g、拉伸应变势能V_e与横向应变势能V_t,经过推导其表达式为

V_g=-

$\int_{0}^{l} [μ + δ (z - l) M_{B} g]$ zdz

(4)

V_e=

$\frac{1}{2} \int_{0}^{l}$ EA

${(v' + \frac{1}{2} w'^{2})}^{2}$ dz

(5)

V_t=

$\frac{1}{2} \int_{0}^{l}$ μ(

$\tilde{l}$ -z)w'²dz

(6)

式中:μ=(ρ_p-ρ_o)Ag 代表立管单位长度湿重,ρ_p、ρ_o分别为立管及立管外部流体的密度,A为立管的横截面积,g 为重力加速度;E为弹性模量;v为立管的轴向形变量;为方便计算可设定 $\tilde{l}$ =l+ $\frac{M_{B} g}{μ}$ . 采用莫里森经验公式,海洋流体力所做的虚功变分为

$\begin{aligned} \delta W_1 & ={ }_0^l Q \delta(r+w) \mathrm{d} z= \\ & { }_0^l\left[-m_{\mathrm{a}}(\ddot{r}+\ddot{w})-\zeta_{\mathrm{d}}\left|v_{\text {rel }}\right| v_{\text {rel }}\right] \delta(r+w) \mathrm{d} z\end{aligned}$

(7)

式中:Q为流体阻力;附加质量m_a=ρ_oA_oC_a,ρ_o为外流密度,A_o为外径横截面面积,附加质量系数C_a=1;阻尼系数ζ_d= $\frac{1}{2}$ ρ_od_oC_d,d_o为立管外径,曳力系数 C_d=1.2;v_rel= $\overset{\cdot}{r}$ + $\overset{\cdot}{w}$ -v_c为立管与海流的相对速度,v_c为洋流速度.同时,母船推力f所做的虚功为

δW₂=fδr

(8)

利用变分法可求出船舶横向运动、立管横向运动和轴向形变方程,其中立管横向运动方程为

$\begin{array}{l}{\left[m+m_{\mathrm{a}}+\delta_{\mathrm{D}}(z-l) M_{\mathrm{B}}\right](\ddot{r}+\ddot{w})+} \\\quad \zeta_{\mathrm{d}}\left|v_{\text {rel }}\right| v_{\mathrm{rel}}-\left(P_{\mathrm{e}} w^{\prime}\right)^{\prime}+ \\\quad E A\left[v^{\prime \prime} w^{\prime}+\left(v^{\prime}+\frac{3}{2} w^{\prime 2}\right) w^{\prime \prime}\right]=0\end{array}$

(9)

式中:δ_D 为狄拉克函数;P_e=μ(l-z)+M_Bg.本研究只考虑横向运动,因此含v(z, t)项舍去.定义一阶小量w'~σ,得:w'²~σ², 只保留一阶小量情况下,横向运动方程最后简化为

$\begin{array}{l}\left(m+m_{\mathrm{a}}\right)(\ddot{r}+\ddot{w})+\zeta_{\mathrm{d}}\left|v_{\text {rel }}\right| v_{\text {rel }}+ \\\mu\left[(z-\widetilde{l}) w^{\prime \prime}+w^{\prime}\right]+ \\\quad \delta_{\mathrm{D}}(z-l) M_{\mathrm{B}}(\ddot{r}+\ddot{w})=0\end{array}$

(10)

基于伽辽金法,可假设

$w(z, t)=\varphi_{i=1}(z) q_{i}(t)=\boldsymbol{\varphi}^{\mathrm{T}}(z) \boldsymbol{q}(t)$

(11)

式中:N为模态函数的数量;q_i(t)为时间坐标且定义q(t)=[ $\begin{array}{l} q_{1} (t) & \dots & q_{N} (t) \end{array}$ ]^T. 模态函数φ_i(z)可以采用相同边界条件下缆索系统(顶端固定,底端自由且附带一个点质量)的固有振型函数且定义φ^T(z)=[ $\begin{array}{l} φ_{1} (z) & \dots & φ_{N} (z) \end{array}$ ].最后得到悬挂立管横向结构动力方程为

$\begin{array}{l}\left[\begin{array}{ll}\boldsymbol{M}_{1} & \boldsymbol{M}_{2}\end{array}\right]\left|\begin{array}{l}\ddot{r} \\\ddot{\boldsymbol{q}}\end{array}\right|+\left[\begin{array}{ll}\boldsymbol{D}_{1} & \boldsymbol{D}_{2}\end{array}\right]\left[\begin{array}{l}\dot{r} \\\dot{\boldsymbol{q}}\end{array}\right]+\boldsymbol{K} \boldsymbol{q}=\boldsymbol{f}_{\mathrm{c}}\\\boldsymbol{M}_{1}={ }_{0}^{l}\left(m+m_{\mathrm{a}}\right) \boldsymbol{\varphi}(z) \mathrm{d} z+M_{\mathrm{B}} \boldsymbol{\varphi}(l) \\\boldsymbol{M}_{2}={ }_{0}^{l}\left(m+m_{\mathrm{a}}\right) \boldsymbol{\varphi}(z) \boldsymbol{\varphi}^{\mathrm{T}}(z) \mathrm{d} z+M_{\mathrm{B}} \boldsymbol{\varphi}(l) \boldsymbol{\varphi}^{\mathrm{T}}(l) \\\boldsymbol{D}_{1}={ }_{0}^{l} \zeta_{\mathrm{d}}\left|v_{\mathrm{rel}}\right| \boldsymbol{\varphi}(z) \mathrm{d} z \\\boldsymbol{D}_{2}={ }_{0}^{l} \zeta_{\mathrm{d}}\left|v_{\mathrm{rel}}\right| \boldsymbol{\varphi}(z) \boldsymbol{\varphi}^{\mathrm{T}}(z) \mathrm{d} z \\\boldsymbol{K}=-\mu_{0}^{l} z \boldsymbol{\varphi}^{\prime}(z) \boldsymbol{\varphi}^{\prime \mathrm{T}}(z) \mathrm{d} z+\mu \tilde{l}_{0}^{l} \boldsymbol{\varphi}^{\prime}(z) \boldsymbol{\varphi}^{\prime \mathrm{T}}(z) \mathrm{d} z \\\boldsymbol{f}_{\mathrm{c}}={0}^{l} \boldsymbol{\varphi}(z) \zeta_{\mathrm{d}}\left|v_{\mathrm{rel}}\right| v_{\mathrm{c}} \mathrm{d} z\end{array}$

(12)

式中:M₁和M₂为质量矩阵;D₁和D₂为阻尼矩阵;K为刚度矩阵;f_c为流体力.

引入日本九州大学Kajiwara实验模型参数^[7]进行仿真计算,预测的立管底端防喷器位置y_d响应曲线如图2所示,与试验数据拟合较好,从而完成模型验证.

图2

图2 立管底端防喷器仿真结果与试验数据对比

Fig.2 Comparison of predicted time trace of LMRP displacement and experimental data

2 控制系统设计

图3所示为基于MPC的立管再入井控制系统原理图.MPC控制器在接收到实船以及立管底端位置信息后,求解下一时刻的最优船速u作为控制指令输入至动力定位系统,从而控制母船及立管的运动轨迹,并记录实船及立管的底端位置作为新的MPC输入.同时在系统中引入观测器对立管受到的洋流扰动进行逼近,以此来提升立管动态响应预测的精度.

图3

图3 基于 MPC的海洋立管再入井作业控制系统原理图

Fig.3 Proposed control system of marine drilling riser in reentry based on MPC

2.1 立管动态响应预测模型

首先将式(12)化为状态空间形式:

$\begin{aligned}\ddot{r}+\boldsymbol{M}_{1}^{-1} \boldsymbol{M}_{2} \ddot{\boldsymbol{q}}+\boldsymbol{M}_{1}^{-1} \boldsymbol{D}_{1} \dot{\boldsymbol{r}}+ \\\quad \boldsymbol{M}_{1}^{-1} \boldsymbol{D}_{2} \dot{\boldsymbol{q}}+\boldsymbol{M}_{1}^{-1} \boldsymbol{K} \boldsymbol{q}=\boldsymbol{M}_{1}^{-1} \boldsymbol{f}_{\mathrm{c}}\end{aligned}$

(13)

设ξ=r+ $M_{1}^{- 1}$ M₂q得:

$\begin{aligned}\ddot{\boldsymbol{\xi}}+ & \boldsymbol{M}_{1}^{-1} \boldsymbol{D}_{2} \boldsymbol{M}_{2}^{-1} \boldsymbol{M}_{1} \dot{\boldsymbol{\xi}}+\boldsymbol{M}_{1}^{-1} \boldsymbol{K} \boldsymbol{M}_{2}^{-1} \boldsymbol{M}_{1} \boldsymbol{\xi}= \\& \boldsymbol{M}_{1}^{-1} \boldsymbol{f}_{\mathrm{c}}+\boldsymbol{M}_{1}^{-1} \boldsymbol{K} \boldsymbol{M}_{2}^{-1} \boldsymbol{M}_{1} r+ \\& \left(\boldsymbol{M}_{1}^{-1} \boldsymbol{D}_{2} \boldsymbol{M}_{2}^{-1} \boldsymbol{M}_{1}-\boldsymbol{M}_{1}^{-1} \boldsymbol{D}_{1}\right) \dot{r}\end{aligned}$

(14)

状态变量设定为X= ${[\begin{array}{l} r & \overset{\cdot}{ξ} & ξ \end{array}]}^{T}$ ,将母船的位移r和立管底端总成的位置y_B作为输出量Y= $[\begin{array}{l} r \\ y_{B} \end{array}]$ ,将母船的速度 $\overset{\cdot}{r}$ 作为控制输入u,得到系统状态空间方程和输出方程为

$\overset{\cdot}{X}$ =

$[\begin{array}{l} 0 & 0 & 0 \\ A_{1} & A_{2} & A_{3} \\ 0 & 1 & 0 \end{array}]$ X+

$[\begin{array}{l} 1 \\ B_{1} \\ 0 \end{array}]$ u+

$[\begin{array}{l} 0 \\ M_{1}^{- 1} f_{c} \\ 0 \end{array}]$

(15)

$\begin{aligned}\boldsymbol{Y}= & \boldsymbol{C} \boldsymbol{X}= \\& {\left[\begin{array}{ccc}1 & 0 & 0 \\1-\boldsymbol{\varphi}(l) \boldsymbol{M}_{2}^{-1} \boldsymbol{M}_{1} & 0 & \boldsymbol{\varphi}(l) \boldsymbol{M}_{2}^{-1} \boldsymbol{M}_{1}\end{array}\right] \boldsymbol{X} }\end{aligned}$

(16)

式中:A₁= $M_{1}^{- 1}$ K $M_{2}^{- 1}$ M₁;A₂=- $M_{1}^{- 1}$ D₂ $M_{2}^{- 1}$ ;A₃=- $M_{1}^{- 1}$ K $M_{2}^{- 1}$ M₁;B₁= $M_{1}^{- 1}$ D₂ $M_{2}^{- 1}$ M₁- $M_{1}^{- 1}$ D₁;C为输出矩阵.设定步长为T,任意第k时间步,已知X(k)和u(k)作为初值,利用欧拉法可对未来H步即k+1~k+H步(H为设定的预测时域步数)状态变量X进行预报,

$\begin{array}{l} X (k + 1) = X (k) + \overset{\cdot}{T X} (k) \\ ︙ \\ X (k + H) = X (k) + T \overset{k + H - 1}{\sum_{i = k}} \overset{\cdot}{X} (i) \end{array}\}$

(17)

2.2 优化函数

在欧拉坐标系Oxz下,设井口坐标(y_d, l),且悬挂立管初始时刻未发生形变时r(0)=0,y(z, 0)=0.设定控制输出目标为Y_d= $[\begin{array}{l} y_{d} \\ y_{d} \end{array}]$ , 表示再入井作业是保证LMRP到达钻井口附近同时需尽量使母船停留在钻井口正上方左右(设置权重矩阵).在任意第k步构建优化函数

$\begin{array}{l} \min \sum_{i=1}^{H}\left(\boldsymbol{Y}(k+i)-\boldsymbol{Y}_{\mathrm{d}}\right)^{\mathrm{T}} \boldsymbol{Q}\left(\boldsymbol{Y}(k+i)-\boldsymbol{Y}_{\mathrm{d}}\right)+ \\\sum_{i=0}^{H-1} R \Delta u(k+i)^{2} \\\text { s. t. } u \in\left[u_{\min }, u_{\max }\right], \quad \Delta u \in\left[\Delta u_{\min }, \Delta u_{\max }\right]\end{array}$

(18)

式中:Q、R 为设定权重;引入输入增量Δu的二次型寻求母船及LMRP到达指定位置,同时母船在再入井作业过程中尽量保持较小速度增量.引入输入指令u及其增量Δu的约束u_min、u_max、Δu_min、Δu_max将母船速度和加速度限制在一定范围内.由于立管预测模型中含有洋流力非线性项,需采用序列二次规划方法对非线性规划问题进行求解,在线求得预测时域内控制输入序列.采用的MPC优化函数结构符合稳定性标准,稳定性证明如下.

定理考虑一般性被控对象离散模型:

$X_{\mathrm{a}}(k+1)=\boldsymbol{f}\left(X_{\mathrm{a}}(k), u_{\mathrm{a}}(k)\right)$

(19)

控制系统在每个周期H内求解的优化命题为

$\begin{array}{l}V_{\mathrm{a}}(k)= \\\quad \min \sum_{i=1}^{H} l_{\text {cost }}\left(X_{\mathrm{a}}(k+i), u_{\mathrm{a}}(k+i-1)\right) \\\text { s. t. } X_{\mathrm{a}} \in X_{\lim }, \quad u_{\mathrm{a}} \in U_{\lim }\end{array}$

(20)

式中:l_cost(X_a, u_a)≥0,且当且仅当X_a=0,u_a=0时l_cost=0.U_lim和X_lim均为包含原点的非空集合.引入终端约束,

X_a(k+H)=0

(21)

假设X_a=0,u_a=0为系统的一个平衡点,并假设每个周期的优化问题都有可行解且能求解得到全局最优解,可认定系统在X_a=0,u_a=0处稳定.

证明利用Lyapunov稳定性进行证明,将V_a(k)作为Lyapunov函数.由于l_cost(X_a, u_a)≥0,可得V_a(k)≥0,仅需证明V_a(k)≥V_a(k+1),假设模型是无偏的,即预测状态与实际状态一致,得:

$\begin{array}{l}V_{\mathrm{a}}(k+1)= \\\quad \min \sum_{i=1}^{H} l_{\text {cost }}\left(X_{\mathrm{a}}(k+i+1), u_{\mathrm{a}}(k+i)\right)= \\\quad \min \left(\sum_{i=1}^{H} l_{\text {cost }}\left(X_{\mathrm{a}}(k+i), u_{\mathrm{a}}(k+i+1)\right)-\right. \\\left.\quad l_{\text {cost }}\left(X_{\mathrm{a}}(k+1), u_{\mathrm{a}}(k)\right)\right)+ \\\quad \min l_{\text {cost }}\left(X_{\mathrm{a}}(k+1+H), u_{\mathrm{a}}(k+H)\right) \leqslant \\\quad-l_{\text {cost }}\left(X_{\mathrm{a}}(k+1), u_{\mathrm{a}}(k)\right)+V_{\mathrm{a}}(k)+ \\\quad \min l_{\text {cost }}\left(X_{\mathrm{a}}(k+1+H), u_{\mathrm{a}}(k+H)\right)\end{array}$

(22)

由于-l_cost(X(k+1), u(k))≤0,且引入终端约束,min l_cost(X(k+1+H), u(k+H))=0.可得V_a(k)≥V_a(k+1),证毕.

2.3 扰动观测器

海洋流体力中引入未知扰动d实时修正,修正后的f_c为

f_c=f_c0+d

(23)

这里定义f_c0= $\int_{0}^{l}$ ζ_d $|v_{r e l}|$ φ(z)v_cdz.设计一款非线性扰动观测器测得 $\hat{d}$ 对真实扰动d进行逼近,设误差e=d- $\hat{d}$ ,假设真实扰动的导数 $\overset{\cdot}{d}$ =0_3×1,则 $\overset{\cdot}{e}$ =- $\overset{\cdot}{\hat{d}}$ .令

$\overset{\cdot}{e}$ =-

$\overset{\cdot}{\hat{d}}$ =-L

$M_{1}^{- 1}$ (d-

$\hat{d}$ )

(24)

L为设计的观测器函数,为使误差e收敛到0,需保证L $M_{1}^{- 1}$ >0,这里令L=( $M_{1}^{- 1}$ )^T,得

$\begin{aligned} \dot{\hat{\boldsymbol{d}}}= & \boldsymbol{L}|| \begin{array}{ccc}0 & 0 & 0 \\ \boldsymbol{A}_1 & \boldsymbol{A}_2 & \boldsymbol{A}_3 \\ 0 & 1 & 0\end{array}|\boldsymbol{X}+| \begin{array}{c}1 \\ \boldsymbol{B}_1 \\ 0\end{array}|u+| \begin{array}{c}0 \\ \boldsymbol{M}_1^{-1} \boldsymbol{f}_{\mathrm{c} 0} \\ 0\end{array}||- \\ & \boldsymbol{L} \boldsymbol{M}_1^{-1} \hat{\boldsymbol{d}}\end{aligned}$

(25)

由于状态变量 $\overset{\cdot}{X}$ 难以测量,令p(X)= $\int_{0}^{t} \overset{\cdot}{L X}$ dt=LX, h= $\hat{d}$ -p(X)得

$\begin{aligned}\dot{\boldsymbol{h}}= & -\boldsymbol{L}|| \begin{array}{ccc}0 & 0 & 0 \\\boldsymbol{A}_{1} & \boldsymbol{A}_{2} & \boldsymbol{A}_{3} \\0 & 1 & 0\end{array}|\boldsymbol{X}+| \begin{array}{c}1 \\\boldsymbol{B}_{1} \\0\end{array}|u+| \begin{array}{c}0 \\\boldsymbol{M}_{1}^{-1} \boldsymbol{f}_{\mathrm{c} 0} \\0\end{array}||- \\& \boldsymbol{L} \boldsymbol{M}_{1}^{-1}(\boldsymbol{h}+p(\boldsymbol{X})) \\\hat{\boldsymbol{d}}= & \boldsymbol{h}+p(\boldsymbol{X})\end{aligned}$

(26)

$\hat{d}$ =h+p(X)

(27)

求解式(26)和(27)可得 $\hat{d}$ ,实现对未知扰动d的实时预测.

2.4 动力定位系统及推进器

根据文献[22],母船采用Abkowitz低频运动模型,在船舶较低航速情况下可对船舶受到的水动力进行简化.在再入井过程中由于船速较低(本研究设定船速范围[-1, 1] m/s),可以忽略二阶速度项.在不考虑船舶艏向角影响下,母船一维横向低频运动可近似表示为

(M-

$M_{\overset{\cdot}{r}}$ )

$\ddot{r}$ -M_r(

$\overset{\cdot}{r}$ -v_c0)=τ+ω

(28)

式中:M为母船的质量; $M_{\overset{\cdot}{r}}$ 为船舶因运动加速度而引起的附加质量;M_r为船舶水动力位置导数;v_c0为水平面洋流流速;τ为推力;ω为零均值白噪声.令u= $\overset{\cdot}{r}$ -v_c0,可构建动力定位系统的控制方程:

$\begin{aligned}\dot{u}= & -\left(M-M_{\dot{r}}\right)^{-1} M_{r} u+\left(M-M_{\dot{r}}\right)^{-1} \tau+ \\& \left(M-M_{\dot{r}}\right)^{-1} \omega\end{aligned}$

(29)

采取PID控制策略可实现动力定位系统模拟.

钻井船一般采用一对全回转导管推进器或吊舱推进器作为主推,可以提供全方向360°的推力,以及一对仅提供横向推力的侧推器.本研究中推进器仅提供船舶纵荡方向的推力,因此仅考虑了单一推进器提供定向推力.引入一阶延时模型用于模拟推进器的动力特性^[23],

$\overset{\cdot}{τ}$ =-A_thr(τ-τ_n)

(30)

式中:τ_n为DPS发出的推力指令;这里只考虑纵荡方向,因此定义A_thr= $\frac{1}{T_{0}}$ , 其中T₀为等效时间常数.

仿真计算流程如下.

(1) 设任意第k时间步(k=1, 2, …),已知立管底端总成位置及母船位置即Y(kT)= $[\begin{array}{l} r (k T) \\ y_{B} (k T) \end{array}]$ (在k=1时,取Y(T)=0_2×1).

(2) 利用式(14)计算得状态变量 $\hat{X}$ (kT)作为X(kT)迭代计算的初值.

(3) 初始化控制输入指令的序列 u₀= ${[\begin{array}{l} u (k T) & u (k T + T) & \dots & u (k T + H T - T) \end{array}]}^{T}$ =0_H_×1.

(4) 利用式(14)可求得状态变量X预估值 $\hat{X}$ = ${[\begin{array}{l} \hat{X} (k T) & \hat{X} (k T + T) & \dots & \hat{X} (k T + H T - T) \end{array}]}^{T}$ . 代入式(15)可求得输出变量Y预估值 $\hat{Y}$ = ${[\begin{array}{l} \hat{Y} (k T) & \hat{Y} (k T + T) & \dots & \hat{Y} (k T + H T - T) \end{array}]}^{T}$ .

(5) 将u₀和 $\hat{Y}$ 代入优化函数即式(17),通过优化方法求得修正的控制序列u₀.

(6) 重复(4)和(5),最终求得最优控制序列 $u_{0}^{*}$ . 将 $u_{0}^{*}$ 中的u(kT)发送给DPS系统并操纵母船并利用立管横向响应模型可求得第k+1步Y(kT+T),然后转入(1)开展下一时间步计算.

3 仿真分析

选用某深海钻井立管系统,立管总长 1 500 m,立管内径和外径分别为0.44 m和 0.5 m,立管密度为 7 900 kg/m³,底端总成净重为100 t.母船净重为 4 200 t,船舶水动力系数 $M_{\overset{\cdot}{r}}$ =-630 000,M_r=-1 260 000.引入东海地区六级海况,三一平均波高 ${\bar{ζ}}_{ω / 3}$ =5.7 m的波浪谱,利用中国船舶及海洋工程设计研究院某船模实验数据响应幅值算子,卷积得到母船纵荡响应(r)如图4所示.悬挂立管受到母船纵荡和v_c=1 m/s均匀洋流作用下的模态响应时间坐标曲线如图5所示,发现立管的第1阶模态始终占据主导,第2阶已被激发,3阶及以上模态分量可以忽略,因此选用模态数N=3是足够的.

图4

图4 波浪谱下船舶纵荡响应

Fig.4 Vessel surge motion in wave spectrum

图5

图5 悬挂立管前4阶模态时间坐标响应曲线

Fig.5 Time traces of first four modal coordinates

设定仿真时间步长T=5 s,预测时域步数H=20,并保持控制时域与预测时域相同.设推进器2 s延迟即A_thr=3 s^-1,权重Q= $[\begin{array}{l} 1 & 0 \\ 0 & 10 \end{array}]$ 和权重R=1,指令约束为-1 m/s<u<1 m/s和-5 m/s<Δu<5 m/s,y_d=50 m即井口坐标(50, 1 500) m,外部均匀洋流且流速为1.0 m/s.采用MATLAB自编程与传统PID再入井系统进行比较.图6为母船与LMRP位置在两种再入井系统作用下的位移(y)响应曲线.PID控制作用下母船和LMRP几乎同时在t≈200 s时越过井口位置进而产生超调,在t≈2 000 s 时,LMRP到达井口位置并保持稳定,从而完成再入井作业.MPC控制作用下,母船能够快速做出响应,悬挂立管也随之快速运动,尽管母船出现较大超调量且在井口上方发生往返运动,在t≈150 s时LMRP到达井口位置并保持稳定.需说明的是,LMRP在t≈120 s 出现波动,但其超调量一直较小,满足工程要求.图7为母船与LMRP位置在两种再入井控制系统作用下的速度( $\overset{\cdot}{y}$ )响应曲线.在PID控制下,LMRP初始阶段受到洋流力作用产生一定波动.MPC控制下,母船快速响应且在t=0~200 s 出现明显波动,因MPC有母船速度和加速度限制,母船速度始终维持在-1~1 m/s.LMRP相较于母船具有一定延迟,但也能随母船快速动作,且在t≈200 s 实现了速度稳定.

图6

图6 母船及立管底端总成位置响应曲线

Fig.6 Time traces of mother vessel and LMRP

图7

图7 母船及立管底端总成速度响应曲线

Fig.7 Time traces of mother vessel velocity and LMRP velocity

再入井作业过程中悬挂立管的第一阶模态通常占主导,一般通过立管的顶端和LMRP横向偏移来探讨悬挂立管弯曲程度.图8为两种再入井系统作用下,立管顶端与LMRP横向偏移时间曲线.PID控制下,立管顶端和LMRP水平距离始终保持较小值,最高不超过8 m. MPC控制下,在t≈60 s时立管顶端和底端距离最高达到40 m,这是由于母船的快速响应和立管因长径比增加致使柔性显著增强造成的.图9为悬挂立管在再入井过程形变响应曲线.图中:l_w为立管与井口的垂直距离.PID控制下,悬挂立管始终保持较小形变.MPC控制下,悬挂立管展现丰富动态响应特征.为深入分析,悬挂立管前两阶模态时间坐标响应如图10所示.PID控制作用下,第1阶模态始终占主导作用.在t=122 s时第1阶和第2阶模态时间坐标同时达到了最高值,再入井过程中悬挂立管形变量较小.在MPC控制作用下, 第2阶模态影响明显增加,在母船和立管回调过程中可能占主导(见t=90 s).这是因为当母船越过钻井口后,需通过DPS对母船实施反方向操纵,而立管由于运动延迟尚未到达井口,激发了悬挂立管第2阶模态.

图8

图8 立管顶端与底端总成之间的横向偏移

Fig.8 Offset of the top-end to the LMRP of riser

图9

图9 再入井过程悬挂立管形变曲线

Fig.9 Deformation curves of hang-off riser in reentry

图10

图10 悬挂立管前两阶模态时间坐标响应曲线

Fig.10 Time traces of first two modal coordinates

分析MPC立管再入系统对洋流力附加模型误差的补偿作用.首先利用wgn函数产生随机高斯白噪声序列,乘以5作为f_c的扰动d,如式(23)所示.施加的扰动和非线性扰动观测器测定的扰动数据如图11所示.受初值及收敛速度影响,t≈150 s 后预测数值与实际扰动数据具有较好贴合度.为防止收敛前的扰动值不能准确预测,可提前启动观测器,然后进行仿真计算.母船与LMRP位置和速度响应曲线分别如图12和图13所示.在扰动影响下,母船和LMRP的响应过程有一定振荡,但依然能很快地完成再入井作业 (LMRP 在t≈120 s到达井口并保持稳定),说明非线性扰动观测器对于扰动具有较好补偿作用.

图11

图11 洋流力附加扰动时间历程曲线

Fig.11 Time trace of added disturbances in hydrodynamic force term

图12

图12 洋流力扰动下母船及LMRP位置响应曲线

Fig.12 Time traces of mother vessel and LMRP under the current resistance disturbance

图13

图13 洋流力扰动下母船及LMRP速度响应曲线

Fig.13 Velocities of mother vessel and LMRP under the current resistance disturbance

分析MPC立管再入系统对洋流速度扰动的补偿作用.利用wgn函数产生0.1 m/s 左右的随机高斯白噪声序列作为均匀流速扰动d_c,即v_c=1.0 m/s+d_c,洋流速度时间曲线如图14所示.母船与LMRP位置和速度的响应曲线分别如图15和图16 所示.扰动作用下,母船和立管的响应过程有一定振荡,但依然能很快地完成再入井作业.底端总成到底钻井口目标位置后,母船仍有一定速度对扰动流速进行动态补偿,保证LMRP始终维持较低速度,有益于安全稳定实现立管再入井作业,从而证明MPC作用下的立管再入系统对于洋流速度扰动具有良好鲁棒性.需指出的是,洋流流速变化导致预测模型中相对速度等参数发生变化,观测器对模型中参数摄动补偿效果有限,仍然需要选择温和海况进行再入井作业.

图14

图14 洋流速度时间历程曲线

Fig.14 Time trace of incoming uniform current

图15

图15 流速扰动下母船及LMRP位置响应曲线

Fig.15 Time traces of mother vessel and LMRP under the flow rate disturbance

图16

图16 流速扰动下母船及LMRP速度响应曲线

Fig.16 Velocities of mother vessel and LMRP under the flow rate disturbance

4 结语

当前海洋钻井工程逐步迈入深海区,恶劣的海洋气候和复杂的海洋环境使钻井立管系统发生脱离(立管底端总成和防喷器断开)的几率大幅增加.如果天气预报成功,可以有时间回收立管并驶离.若不能及时预报则必须进行紧急脱离.天气转好后,则需要将底端总成和防喷器重新连接,称为再入井作业.由于深海区天气与海况复杂多变,确认海洋环境适合时需要尽快完成再入井作业.但是悬挂立管系统因长径比的大幅增加,立管的柔性显著增强,这对在母船激励和复杂海况下(即使海况温和)安全快速完成再入井作业提出新的严峻挑战.

基于修正哈密顿原理建立底端含集中质量(底端总成)的柔性悬挂立管系统仿真模型,结合立管动态响应预测模型及井口位置设计优化函数和约束,构建非线性扰动观测器实现对洋流力的模型误差和洋流速度的扰动补偿,尝试提出一种基于MPC的深海钻井立管再入井控制系统.仿真发现:相比于传统钻井立管再入井作业PID控制,在MPC控制系统作用下,母船和悬挂立管可快速做出响应,立管系统能够安全稳定地实现再入井作业,能较好地处理洋流力模型误差问题,且在洋流速度扰动下具有良好鲁棒性.

本研究属于寻找海洋立管快速再入井策略初步探索,有许多不足之处.例如,只考虑了立管横向(母船纵荡)运动;当洋流泄涡频率接近系统固有频率时,立管可能发生涡激振动.下一步将开展横向和纵向耦合动态响应分析及二维再入控制系统研究,此时需要分析推进器布置及推力分配问题.同时为防止立管形变较大产生破坏,可在优化函数中引入悬挂立管最大形变约束,这些将在后续工作中逐步完成.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

MURAI

, YAMAMOTO

An experimental analysis of the internal flow effect on marine risers

[C]// Proceedings of the Internal Conference on Marine Technology. Dhaka, Bangladesh: MARTEC, 2010: 159-165.

[2]

王宗义, 肖坤, 庞永杰,

等.

船舶动力定位的数学模型和滤波方法

[J]. 哈尔滨工程大学学报, 2002, (4): 24-28.

WANG

Zongyi

, XIAO

Kun

, PANG

Yongjie

et al.

Mathematical models and filtering methods for dynamic ship positioning

[J]. Journal of Harbin Engineering University, 2002(4): 24-28.

[3]

SUZUKI

, YOSHIDA

, SHIDA

et al.

Active control of riser deformation and vessel motion for automatic entry/reentry system

[C]// Proceedings of 4th International Offshore and Polar Engineering Conference. Osaka, Japan: ISOPE, 1994: 216-223.

[4]

SUZUKI

, TAO

, WATANABE

Learning tracking controllers under unknown disturbances for the installation of rigid and flexible structures

[J]. Journal of Marine Science & Technology, 1999, 4(4): 187-199.

[5]

TAKAYUKI

, KAZUHISA

, KAJIWARA

et al.

On vibration control of flexible pipe in ocean drilling system

[C]// Proceedings of the Sixteenth International Offshore and Polar Engineering Conference. San Francisco, USA: ISOPE, 2006: 26-30.

[6]

KOTERAYAMA

, NAKAMURA

, KAJIWARA

et al.

Development of a training simulator for dynamic reentry operations of a riser pipe hanged off

[C]// Proceedings of the Eighteenth International Offshore and Polar Engineering Conference. Vancouver, Canada: ISOPE, 2008: 67-70.

[7]

KOTERAYAMA

, KARJIWARA

, NAKAMURA

et al.

Model experiments of reentry control based on riser pipe dynamics

[C]// The Eighth ISOPE Pacific/Asia Offshore Mechanics Symposium. Bangkok, Thailand: ISOPE, 2008: 8-12.

[本文引用: 5]

[8]

KARJIWARA

, NORIDOMI

Reentry control system design for a riser pipe experimental model under steady current

[C]2009 ICCAS-SICE. Fukuoka, Japan: IEEE, 2009: 3896-3901.

DOI:10.1007/s13344-012-0048-7 URL [本文引用: 3]

[9]

Y L

, YAO

B H

, ZHENG

et al.

Research on marine riser in different installations stages of subsea production tree

[C]// Shanghai, China. China: IEEE, 2016: 920-922.

[本文引用: 2]

[10]

WANG

S W

, XU

X S

, YAO

B H

et al.

A finite difference approximation for dynamic calculation of vertical free hanging slender riser in re-entry application

[J]. China Ocean Engineering, 2012, 26(4): 637-652.

[11]

王盛炜, 徐雪松, 姚宝恒,

等.

针对立管重入井过程中等路径段加速度的改进型蚁群优化算法

[J]. 上海交通大学学报, 2013, 47(10): 1585-1590.

[本文引用: 3]

WANG

Shengwei

, XU

Xuesong

, YAO

Baoheng

et al.

Improved ant colony optimization for accelerations of equal-length paths of riser re-entry

[J]. Journal of Shanghai Jiao Tong University, 2013, 47(10): 1585-1590.

[本文引用: 3]

[12]

WANG

S W

, XU

X S

, LU

Movement optimization of freely-hanging deepwater risers in reentry

[J]. Ocean Engineering, 2016, 116: 32-41.

DOI:10.1016/j.oceaneng.2016.02.029 URL [本文引用: 3]

[13]

TRIANTAFYLLOU

M S

Cable mechanics for moored floating systems

[C]// Proceedings of the Seventh International Conference on the Behavior of Offshore Structures. Massachusetts, USA: Pergamon, 1994: 57-79.

DOI:10.1016/j.oceaneng.2008.01.006 URL [本文引用: 1]

[14]

CHATJIGEORGIOU

I K

A finite differences formulation for the linear and nonlinear dynamics of 2D catenary risers

[J]. Ocean Engineering, 2008, 35(7): 616-636.

[15]

Y L

, CAO

J J

, YAO

B H

et al.

Dynamic behaviors of a marine riser with variable length during the installation of a subsea production tree

[J]. Journal of Marine Science and Technology, 2017, 23(3): 1-11.

DOI:10.1007/s00773-017-0462-x URL [本文引用: 1]

[16]

胡永利.

深海钻井立管重入井作业运动特性及其优化研究[D]. 上海: 上海交通大学, 2018.

Yongli

Study of properties and optimization of deep-sea drilling riser in reentry[D]. Shanghai: Shanghai Jiao Tong University, 2018.

[17]

Y L

, YAO

B H

, LIAN

Analysis of drilling riser in re-entry operation under complex currents

[J]. Journal of Ship Mechanics, 2019, 23(12): 1528-1540.

[18]

YAMAMOTO

, MOROOKA

C K

, UENO

Dynamic behavior of a semi-submersible platform coupled with drilling riser during re-entry operation in ultra-deep water

[C]// Proceedings of the ASME 2007 26th International Conference on Offshore Mechanics and Arctic Engineering. San Diego, California, USA: ASME, 2007: 239-248.

[19]

YAMAMOTO

An automatic approaching control for the petroleum well re-entry operation in ultra-deep water

[C]// Proceedings of the 8th IFAC International Conference on Maneuvering and Control of Marine Craft. Guarujá, Brazil: IFAC, 2009: 388-393,

DOI:10.1016/j.apor.2018.11.013 [本文引用: 1]

[20]

YAMAMOTO

, MOROOKA

C K

Feedback control system for blow-out preventer positioning

[J]. Applied Ocean Research, 2019, 82: 362-369.

During the drilling of ultra-deep-water subsea petroleum wells, a blow-out preventer (BOP), a piece of safety equipment, must be assembled on the wellhead. The BOP is suspended using the drilling riser during the wellhead approach operation, and the riser's top end is connected to the floating platform rig. This article presents a feedback control system for the automatic approach of the BOP to the wellhead. Compared to state-of-the-art controls, ours does not require ancillary thrusters installed alongside the riser nor inclination sensors atop of the drilling riser. Additionally, our proposed control embeds a closed-loop dynamic positioning system, thus retaining the characteristics of the original control system and adding an extra closed-loop. This eases implementation of the BOP approach control to an existing platform. To calculate the optimal gains for the BOP controller, we assume a linear system for the riser, including only the pendulum-shape. The simulation is carried out using nonlinear models for both riser and floating platform. We assume an International Towing Tank Conference standard semi-submersible platform, coupled with a 3000-m free-hanging vertical riser for the time-domain simulation. The results show the BOP tracking to be a step-shaped input signal under current and wave loads. A discussion of the performance of feedback control under different environmental loads is also included.

[21]

毛良杰, 付燈煌, 曾松.

正常连接与悬挂撤离工况下深水钻井隔水管动力特性与安全分析

[J]. 船舶力学, 2022, 26(2): 250-263.

MAO

Liangjie

, FU

Denghuang

, ZENG

Song

Dynamic characteristics and safety analysis of deep water drilling riser under normal connected condition and hang-off evacuation condition

[J]. Journal of Ship Mechanics, 2022, 26(2): 250-263.

[22]

FOSSEN

T I

Guidance and control of ocean vehicles

[M]. Chichester, England: John Wiley & Sons, 1994.