预条件GMRES(m)算法在钻井隔水管力学分析中的应用

李朝玮,许亮斌,盛磊祥,李梦博,姜智博

(中海油研究总院有限责任公司 技术研发中心,北京 100028)

在深水浮式钻井中,钻井隔水管(钻井立管)系统顶端与钻井平台/船相接[1],底端与水下井口连接(图1),在海水载荷激励下和浮式平台运动的牵引下,隔水管容易发生侧向偏移、弯曲变形以及水平和纵向振动．如果隔水管受力变形和振动幅度超出限度,将会妨碍钻完井作业,诱发隔水管弱点位置疲劳破坏,影响整个钻井系统的稳定和安全．因此,在深水钻井设计时,需要对钻井隔水管系统进行详细的力学分析．工程设计中对钻井隔水管力学分析一般依靠ANSYS、ABAQUS、OrcaFlex和DeepRiser等国外软件,国内开发的隔水管分析软件尚未达到成熟的商业应用阶段．为了开发和优化国产化隔水管软件,均需编程求解隔水管力学方程．

虽然隔水管力学问题较多(包括静态/准静态分析[2],波激[3]、涡激[4]、参激[5]振动分析,应急解脱及反冲响应分析[6],悬挂撤离动态分析[7],隔水管-井口-导管耦合力学分析[8],干涉与碰撞分析[9]等),但求解这些问题时,一般均要将隔水管控制方程通过有限元法或有限差分法进行离散,形成一个大型稀疏非对称非正定的线性方程组;求解每个时间步的方程组未知量即各自由度的变化,进一步可得到管体的受力变形和动态响应．

在深水中,更长的隔水管需要更多的离散单元,增加了离散方程组的阶数;考虑管体、外力以及边界条件的非线性因素后,在方程组的矩阵中可能引入超大元素或降低系数矩阵的稀疏性．这些极大地增加了隔水管离散方程组高精度高效率求解的难度．目前,研究者们的关注点多在隔水管动力学方程形式、外激力的求解方法和动态响应特征分析,未见文献系统分析隔水管力学离散方程组的高精度高效率求解问题．但是,采用常规的线性方程组求解方法计算这类大型方程组时,在数据存储、计算效率和求解精度方面均不理想．

文中在对比分析了几种大型线性方程组的求解方法后,优选出了“预条件的重开始广义极小残量法(iLU-GMRES(m),incomplete LU factorization-generalized minimal residual method (with restart))”,以“海流+表面波+内孤立波+平台纵荡运动”共同作用下的1 080 m长隔水管为例,采用Matlab编程进行了GMRES算法的应用和隔水管动力学分析．

1 钻井隔水管动力学计算方法

钻井隔水管力学分析用于计算管体在波流等外载荷激励下关键节点的位移、转角、弯矩、剪力、应力等随时间的变化,确定管体力学参数的极值窗口,并为弱点疲劳分析提供时程数据．

1.1 控制方程的时域动态求解

如图1,设z为纵坐标,垂直向上为正,m;y为管体水平方向位移,顺流向为正,m;t为时间,s．采用文献[7]中的方法建立钻井隔水管动力学控制方程并进行有限元离散,得到任意时刻t的动力学有限元方程:

(1)

图1 深水钻井隔水管系统示意Fig.1 Schematic diagram of a deepwaterdrilling riser system

(2)

其有效刚度矩阵和有效载荷向量分别为:

(3)

1.2 隔水管力学方程求解的计算量分析

图2显示了该方程组的系数矩阵的非零元素分布(为便于图示,此图仅显示系数矩阵两个端部和中间共20个节点的数据),是一个带宽为6的带状稀疏矩阵;中间位置节点的系数矩阵值变化相对平缓,但上下边界附近的矩阵元素值波动剧烈,会增加数值计算收敛过程的难度;由于考虑了有效轴向力和阻尼等因素,并引入了边界条件,破坏了系数矩阵的对称性,对称位置上的部分元素有较小的偏差,所以系数矩阵是非对称非正定的．考虑内孤立波作用时间范围内共4 000 s下的隔水管动力学响应,每个时间步取2 s,则在整个动力学分析中需求解2 000次1 202×1 202的大型稀疏非对称线性方程组．

如果采用三维离散、更长的分析模型、更密集的单元和更小的时间步,则线性方程组的阶数和时间步都会大量增加．在进行钻井隔水管动力学分析时,由于数值计算的误差来源、计算时间的耗费以及中间数据的内存占用量主要来自每个时间步求解一次这一大型稀疏非对称线性方程组,所以选择迭代时间短、计算精度高和内存占用量小的线性方程组求解算法至关重要．

表1 钻井隔水管动态分析实例数据

图2 系数矩阵的非零元素分布Fig.2 Distribution of nonzero elementsin the coefficient matrix

2 大型稀疏非对称线性方程组求解算法优选

2.1 线性方程组求解算法对比分析

线性方程组的解法超过20余种,可分为直接法、线性定常迭代法和非定常迭代法[15-18]．对于大型稀疏线性方程组,直接法(如Gauss消去法、LU分解法等)将使系数矩阵丧失稀疏特性,极大地增加计算机的存储量,并导致无法求出收敛解;线性定常迭代法(如Gauss-Seidel迭代法、超松弛迭代法SOR等)在迭代收敛速度和求解精度方面也不理想．

近半个世纪以来,学者们针对不同类型的大型线性方程组,开发了十余种非定常迭代法,Matlab软件收录了其中的10种成熟算法,如表2．由于隔水管力学分析系数矩阵为大型稀疏非对称非正定的,可选用的方法只有BiCGSTAB(或BiCG、BiCGSTABL)、CGS、GMRES、LSQR、QMR、TFQMR．

表2 Matlab中的线性方程组迭代算法对系数矩阵的适用性分析

图3 6种迭代算法的迭代次数和相对残差Fig.3 Diagram of iteration numbers-relative residual among six iteration methods

2.2 重开始广义极小残量法GMRES(m)

广义极小残量法(GMRES)[17-18]对求解大型稀疏非对称非正定线性方程组有良好的数值稳定性和较高的计算效率,且在科学计算和工程分析中应用较多[19-20]．在文献[17]中,极小化处在Krylov子空间的解的残量范数,通过构造正交向量序列并用其线性组合来更新解向量,其中线性组合的系数通过求解一个较小规模的最小二乘问题而得．GMRES算法每迭代一步都要进行Arnoldi过程．随着子空间维数的增大,其计算量和存储量都会增加．为了克服这一“长拖”问题,改善条件数以加速收敛,Saad和Schultz提出了重开始(restart)的GMRES算法,记为GMRES(m),其中，m为重新开始参数,m取值通常远远小于线性方程组的阶数．图4为GMRES(m)的算法步骤流程,详细内容可参考文献[17-20]．

2.3 预条件处理

对于大型线性方程组Ax=b,通常对其直接迭代求解的收敛速度很慢,这跟系数矩阵A的条件数较大而趋于病态直接相关．为了提高迭代求解的收敛速度,一般需对系数矩阵A进行预条件处理,即使得系数矩阵趋向于单位矩阵,特征值聚到一个很小的复数域之中．取预处理矩阵M,有两种预条件处理方式:

图4 GMRES(m)算法流程Fig.4 Flow chart of GMRES(m) algorithm

常用的预条件技术包括不完全LU分解(iLU)和对角线型不完全LU分解(DiLU)[21]等．文中采用较简单的iLU分解法,即将系数矩阵A分解为一个上三角矩阵和一个单位下三角矩阵乘积的近似形式:A=LU+R,R为满足某些极限条件的残留矩阵．则预处理矩阵取为M=LU．

图5 预条件处理后的迭代次数和相对残差Fig.5 Diagram of iteration numbers-relative residual among six iteration methods after preconditioning

3 方法对比和应用实例

3.1 方法对比

分别采用① ANSYS模拟、② DeepRiser软件模拟、③ 文中程序和LU分解法、④ 文中程序和iLU-GMRES(m)算法,以表1中的数据为基础进行连接工况下的隔水管动态计算,对比计算结果见表3．对比方法③ 中,选择LU分解法是因为LU分解法是线性方程组的较常规解法．由于在商业软件中不易添加内孤立波载荷,方法对比先只考虑海流(OC)、表面波(SW)和钻井平台纵荡运动(VS)的组合,其中平台运动假设钻井平台定位良好,静态漂移量为零,考虑了平台的长期漫漂运动和随不规则波浪的瞬时响应．

从表3对比结果可看出,采用文中程序和iLU-GMRES(m)算法计算的隔水管动态响应的横向偏移最大位移及对应的水深位置与ANSYS、DeepRiser软件模拟的结果基本一致,但采用LU分解法时虽然也能得到结果,但由于计算时间步较多,累积了较大误差,其结果与其他3种方法差别很大．在计算效率上,采用文中程序和iLU-GMRES(m)算法时,由于每个时间步求解线性方程组时只需耗时约0.007 s,动态模拟的总耗时显著优于DeepRiser和ANSYS软件模拟．

表3 方法对比(OC+SW+VS)

3.2 应用实例

进一步考虑内孤立波(IW)与其他载荷的叠加作用,以表1中的数据为基础,采用文中程序和iLU-GMRES(m)算法进行连接工况的隔水管动态计算．其中iLU-GMRES(m)算法的重开始次数m取为6,内孤立波计算采用两层流的KdV模型[2, 3, 7]．

图6为海流、表面波、内孤立波、平台运动共同作用下的钻井隔水管动态响应时程曲线．如图6(a～b),各时刻点的各深度位置动态响应曲线的边界构成了动态响应参数的极值包络线,反映了动态响应的最危险情况．图6(c～f)为隔水管关键位置点的参数时程曲线,图6(d)中,水下100 m(Mises应力的一个极值点位置)处的应力由平衡位置的85 MPa扩大到峰值处的95 MPa;图6(e)中,顶部挠性接头的转动角度幅值从(-0.6°～0.2°)扩大到(-1.0°～0.3°);图6(f)中,底部挠性接头的转动角度幅值从(-0.4°～0.2°)扩大到了(-0.7°～0.5°)．

图7为不同载荷组合时的钻井隔水管性能参数动态响应的极值包络线对比曲线,进一步显示了内孤立波显著增加了钻井隔水管的顺流向偏移、弯矩、剪力、Mises应力的幅值,如图7(c)中内孤立波使管体剪力在波界面深度附近处出现了一个峰值．此外,由于文中考虑了顶底挠性接头的旋转刚度,所以计算的弯矩、剪力和应力均在顶部、底部出现大幅增加．

中国南海内孤立波分布广泛,发生频繁[22]．为了减弱内孤立波对钻井隔水管系统的影响,需做好海上油气作业区的内孤立波监测工作,确定内孤立波的来流方向、流速以及海水分界面的深度范围和海水密度差异．同时,对上层海水范围内以及内孤立波分界面一定深度范围内的隔水管选型,建议选用更大强度的隔水管(增加单根壁厚或钢级)以抵抗内孤立波叠加海流作用的强剪切破坏,并在这些深度范围内尽量选用裸单根以降低管体受到的流体拖曳力．

图6 多种载荷共同作用下的连接状态隔水管性能参数动态响应Fig.6 Dynamic responses of connected riser properties under combined loads

图7 多种载荷共同作用下的连接状态隔水管性能参数动态响应Fig.7 Dynamic responses of connected riser properties under combined loads

4 结论

(1) 钻井隔水管力学计算和分析是深水钻井工程设计的重要方面．在对隔水管结构动力学进行数值求解时,离散后的控制方程将转化为一个大型稀疏非对称线性方程组．文中优选了“预条件处理的重开始广义极小残量法(iLU GMRES(m))”,结合Wilson-θ法进行深水钻井隔水管结构动力学数值求解,单个时间步迭代6次就能获得相对残差接近7.01×10-20的解,单次迭代耗时只需0.007s,减少了迭代时间,降低了中间数据的内存占用量,提高了计算效率．

(2) 采用该套算法对海流、表面波、内孤立波和平台纵荡运动共同作用下的钻井隔水管动态性能进行计算,揭示出内孤立波显著扩大了隔水管的顺流向偏移、转角、弯矩、剪力和应力的极值窗口,增加了隔水管系统的危险性．在南海油气开发时,建议做好作业区的内孤立波监测,并在上层海水和分界面附近一定深度范围内使用隔水管裸单根并选择更大钢级和壁厚的隔水管．